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PRELIMINARY  DEVELOPMENT  OF  A  THREE-DIMENSIONAL  NONLINEAR 
APPROXIMATION  PROCEDURE  FOR  TURBOMACHINERY  DESIGN 
OPTIMIZATION  APPLICATIONS 


Stephen  S.  Stahara.  Jasopin  Lee, 
and  John  R.  Spreiter 


SUMMARY 

A  theoretical  investigation  was  carried  out  involving  the  preliminary 
development  of  a  novel  three-dimensional  nonlinear  approximation  method 
capable  of  rapidly  determining  approximate  solutions  for  highly  nonlinear  flows 
such  as  typically  occur  in  turbomachinery  applications.  The  ultimate  goal  of  the 
study  was  the  preliminary  demonstration  of  the  accuracy  and  utility  of  such  a 
method  in  the  three-dimensional  turbomachinery  optimization  design 
environment.  The  specific  objectives  of  the  present  investigation  were  the 
theoretical  development  of  the  three-dimensional  approximation  procedure  for 
highly  sensitive,  nonlinear  flows  involving  multiple  shock  waves,  combination  of 
the  approximation  method  with  a  nonlinear  3-D  supercritical  transonic  flow 
solver,  coupling  of  the  combined  approximation  method  and  nonlinear  flow 
solver  with  an  optimization  design  procedure,  and  finally  testing  of  the  complete 
approximation/3-D  flow  solver/optimization  code  on  problems  relevant  to  the 
three-dimensional  turbomachinery  design  optimization  environment. 


1.  INTRODUCTION 


The  remarkable  success  of  advanced  computational  methods  for  determining 
accurate  solutions  to  increasingly  complex  fluid  dynamic  phenomena  has  now 
been  well  established  across  a  broad  range  of  flow  problems  of  practical 
interest.  What  has  also  become  simultaneously  apparent  with  this  success  is 
that  a  major  impediment  exists  to  the  implementation  of  these  emerging  codes 
when  employed  in  highly-repetitive  usage  applications.  This  is  due  to  the 
excessive  computational  demands  required  by  their  straightforward  application. 
Many  such  applications  simply  cannot  afford  the  computational  cost  associated 
with  the  repetitive  use  of  these  high-level  numerical  solvers.  Consequently,  a 
need  clearly  exists  for  the  development  and  implementation  of  complimentary 
nonlinear  approximation  methods  that  are  sufficiently  general  and  accurate  to 
be  used  in  conjuction  with  these  advanced  codes  to  reduce  their  computational 
demands.  While  this  need  exists  across  a  spectrum  of  aerodynamic  and  fluid 
dynamic  uses,  it  is  particularly  high  in  turbomachinery  applications.  For  that 
application,  both  the  basic  fluid  dynamic  computation  is  highly  time  consuming 
and  the  number  of  variable  flow  and  geometry  parameters  are  large,  resulting  in 
any  turbomachinery  parametric  or  design  study  being  computationally 
expensive  under  the  best  of  circumstances,  and  in  many  instances  using  more 
advanced  codes  prohibitively  so. 

The  ultimate  objective  beyond  this  initial  study  is  to  develop  and  demonstrate 
the  feasibility  of  such  approximation  methods  for  substantially  reducing  the 
overall  computational  requirements  necessary  for  general  3-D  turbomachinery 
design  or  parametric  optimization.  It  is  conceived  that  these  approximation 
methods  would  be  coupled  with  high  run-time  3-D  turbomachinery 
computational  flow  solvers,  and  would  be  used  in  conjunction  with  these 
solvers  in  applications  where  large  numbers  of  related  solutions  are  needed. 
The  computational  time  saving  would  be  accomplished  by  employing  these 
rapid  approximation  methods  to  decrease  to  a  minimum  the  actual  number  of 
expensive  3-D  flow  solutions  needed  in  any  optimization  or  design  study.  The 
actual  implementation  would  entail  coupling  the  rapid  approximation  method 
with  a  3-D  tubomachinery  flow  solver  and  a  design  optimization  procedure  into 
a  combined  code,  and  then  employing  the  approximation  method  in  the 
combined  code  together  with  a  certain  minimum  number  of  pre-determined  3-D 
flow  solutions  to  then  predict  all  of  the  flow  solutions  subsequently  required  by 
the  optimization  driver  as  that  procedure  searches  through  the  design  variable 
solution  space  to  reach  the  final  optimum  design  point. 

That  such  procedures  are  in  fact  achievable  has  been  successfully 
demonstrated  for  two-dimensional  flows.  In  studies  made  by  the  present  author 
and  reported  in  Refs.  1-7,  a  remarkable  nonlinear  approximation  method  was 
developed  and  extensively  tested  on  a  wide  range  of  both  continuous  and 
discontinuous  nonlinear  flow  problems.  Its  ability  to  accurately  predict 
nonlinear  solutions  of  primary  interest  to  this  study  was  first  confirmed  in  case 
studies  involving  a  variety  of  strongly  nonlinear  transonic  flows.  The  method 
was  then  coupled  with  an  optimization  procedure  and  tested  on  several  two- 
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dimensional  design  problems.  The  results  demonstrated  the  potential  of  the 
approximation  method  to  reduce  the  computational  work  in  such  applications  by 
an  order  in  magnitude  with  no  degradation  in  accuracy  of  the  final  optimized 
design  result. 

The  theoretical  basis  of  the  approximation  procedure  encompasses  both 
inviscid  and  viscous  flows.  2-D  and  3<D  situations,  and  both  steady  and 
unsteady  flows.  It  has  been  developed  most  extensively  for  the  steady  situation, 
however,  and  primarily  focused  in  initial  design  applications  toward 
design/optimization  studies  of  highly-nonlinear  2-D  internal  and  external 
transonic  flows.  A  notable  point  regarding  the  development  and  use  of  the 
approximation  method  for  the  application  described  here  is  that  the  method 
does  not  depend  on  any  particular  turbomachinery  flow  analysis  code. 
Consequently,  obsolescence  of  the  methodology  due  to  future  analysis  code 
development  is  not  a  factor. 

The  work  reported  here  involves  the  initial  development  and  extension  of  these 
methods  and  concepts  to  the  three-dimensional  turbomachinery  optimization 
design  problem.  The  specific  implementation  involves  development  of  the 
nonlinear  approximation  method  in  a  form  suitable  for  predicting  surface 
properties  on  three-dimensional  turbomachinery  blades;  and  then  integration  of 
that  form  of  the  approximation  method  with  a  3-D  design  optimization 
procedure.  The  QNMDIF  optimization  procedure  recently  developed  at 
NASA/Ames  Research  Center  and  successfully  demonstrated  in  a  variety  of 
aerodynamic  design  optimization  applications  (Ref.  8)  was  selected  for 
implementation  in  this  study.  That  procedure  consists  of  the  QNMDIF 
optimization  driver  (Ref.  9)  coupled  with  the  TWING  three-dimensional  full 
potential  solver  (Ref.  10)  for  determining  supercritical  transonic  flows  past  wing 
or  isolated  3-D  blade  geometries. 
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2.  ANALYSIS 


2.1  Description  of  the  Nonlinear  Approximation  Concept 

The  classical  approach  of  developing  a  perturbation  or  approximation  analysis  ■ 
that  is,  by  establishing  and  solving  a  series  of  linear  perturbation  equations  in 
the  manner  elegantly  described  arid  illustrated  by  Van  Dyke  (Ref.  1 1)  ~  appears 
as  an  obvious  choice  for  the  current  application.  However,  results  from  the  work 
reported  in  Ref.  1  demonstrate  that  for  applications  to  sensitive  flows  such  as 
occur  in  most  turbomachinery  situations,  the  basic  linear  variation  assumption 
fundamental  to  such  a  technique  is  sufficiently  restrictive  that  the  permissible 
range  of  parameter  variation  is  so  small  to  be  of  little  practical  use.  An 
interesting  and  novei  alternative  to  the  linear  perturbation  equation  approach 
has  recently  been  successfully  examined  in  which  a  correction  or 
approximation  technique  is  used  that  employs  two  or  more  nonlinear  base 
solutions.  For  the  approximation  method,  the  basic  solution  is  determined 
simply  by  differencing  two  nonlinear  flow  solutions  removed  from  one  another 
by  some  nominal  change  of  a  particular  flow  or  geometrical  quantity.  A  unit 
approximation  solution  is  then  obtained  by  dividing  that  result  by  the  change  in 
the  varied  quantity.  Related  solutions  are  determined  by  multiplying  the  unit 
approximation  solution  by  the  desired  parameter  change  and  adding  that  result 
to  the  base  flow  solution.  This  simple  procedure,  however,  only  works  directly 
for  continuous  flows  for  which  the  parameter  change  does  not  alter  the  solution 
domain.  For  those  parameter  changes  which  change  the  flow  domain, 
coordinate  stretching  is  necessary  to  ensure  proper  definition  of  the  unit 
^proximation  solution.  Similarly,  for  discontinuous  flows,  coordinate  straining 
is  necec^sary  to  account  for  movement  of  discontinuities  due  to  the  parameter 
change.  We  will  discuss  in  detail  the  importance  of  coordinate  straining  to  the 
approximation  method  below. 

The  attractiveness  of  such  an  approximation  method  is  that  it  is  not  restricted  to 
a  linear  variation  range  but  rather  replaces  the  nonlinear  variation  between  two 
base  solutions  with  a  linear  fit.  This  de-emphasizes  the  dependence  and 
sensitivity  inherent  in  the  linear  perturbation  equation  method  on  the  local  rate 
of  change  of  the  base  flow  solution  with  respect  to  the  varied  quantity.  For  many 
applications,  particularly  at  supercritical  transonic  speeds,  the  flow  is  highly 
sensitive,  and  the  linear  range  of  parameter  variation  can  be  sufficiently  small  to 
be  of  no  practical  use.  Furthermore,  other  than  the  approximation  of  a  linear  fit 
between  two  nonlinear  base  solutions,  this  new  method  is  not  restricted  by  any 
further  approximations  with  respect  to  the  governing  differential  equations  and 
boundary  conditions.  Rather,  it  retains  the  full  character  of  the  original  methods 
used  to  calculate  the  base  flow  solutions.  Most  importantly,  no  perturbation 
differential  equations  have  to  be  posed  and  solved,  only  algebraic  ones.  In  fact, 
it  isn't  even  necessary  to  know  the  exact  form  of  the  perturbation  equation,  only 
that  it  can  be  obtained  by  some  systematic  procedure  and  that  the  perturbations 
thus  defined  will  behave  in  some  'generally  appropriate'  fashion  so  as  to  permit 
a  logical  perturbation  analysis.  For  situations  involving  variations  of  physical 
parameters,  such  as  reported  here,  the  governing  perturbation  equations  are 
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usually  transparent,  or  at  least  readily  derivable.  Finally,  in  applying  this 
method  it  isn't  necessary  to  work  with  primitive  variables;  rather  the  procedure 
can  be  applied  directly  to  the  final  solution  quantity  desired.  An  important 
qualification  of  this  method  is  that  the  two  base  solutions  required  for  each 
parameter  change  considered  must  be  topologically  similar,  i.e..  discontinuities 
or  other  characteristic  features  must  be  present  in  both  base  solutions  used  to 
establish  the  unit  approximation  solution. 

The  concept  of  employing  coordinate  straining  to  remove  nonuniformities  from 
perturbation  solutions  of  nonlinear  problems  is  well  established  and  was 
originally  suggested  by  Ughthill  (Ref.  12)  four  decades  ago.  The  basic  idea  of 
the  technique  is  that  a  straightfonward  perturbation  solution  may  possess  the 
appropriate  form,  but  not  quite  at  the  appropriate  location.  The  procedure  is  to 
strain  slightly  the  coordinates  by  expanding  them  as  well  as  the  dependent 
variables  in  an  asymptotic  series.  It  is  often  unnecessary  to  actually  solve  for 
the  straining,  it  generally  can  be  established  by  inspection.  The  final  uniformly 
valid  solution  is  then  found  in  implicit  form,  with  the  strained  coordinate 
appearing  as  a  parameter. 

In  the  original  applications  of  the  method  (Ref.  11),  it  was  applied  in  the 
'classical'  sense;  that  is.  series  expansions  of  the  dependent  and  independent 
variables  in  ascending  powers  in  some  small  parameter  were  inserted  into  the 
full  governing  equation  and  boundary  conditions,  and  the  individual  terms  of  the 
series  determined.  An  ingenious  variation  in  the  application  of  the  method  was 
made  by  Pritulo  (Ref.  13)  who  demonstrated  that  if  a  perturbation  solution  in 
unstrained  coordinates  has  been  determined  and  found  to  be  nonuniform,  the 
coordinate  straining  required  to  render  that  solution  uniformly  valid  can  be 
found  by  employing  straining  directly  in  the  known  nonuniform  solution,  and 
then  solving  algebraic  rather  than  differential  equations.  The  idea  of 
introducing  strained  coordinates  a  posteriori  has  since  been  applied  to  a  variety 
of  different  problems  (Ref.  1 1 )  and  forms  the  basis  of  the  current  application. 

The  fundamental  idea  underlying  coordinate  straining  as  it  relates  to  the 
application  of  perturbation  or  approximation  methods  to  nonlinear  flows  as  we 
apply  them  here  is  illustrated  geometrically  in  Figure  1.  In  the  upper  plot  on  the 
left,  two  typical  transonic  pressure  distributions  are  shown  for  a  highly 
supercritical  flow  about  a  nonlifting  symmetric  profile.  The  distributions  can  be 
regarded  as  related  nonlinear  flow  solutions  separated  by  a  nominal  change  in 
some  geometric  or  flow  parameter.  The  shaded  area  between  the  solutions 
represents  the  perturbation  result  that  would  bo  obtained  by  directly  differencing 
the  two  solutions.  We  observe  that  the  perturbation  so  obtained  is  small 
everywhere  except  in  the  region  between  the  two  shock  waves,  where  it  is  fully 
as  large  as  the  base  solutions  themselves.  This  clearly  invalidates  the 
perturb^ion  technique  in  that  region  and  most  probably  somewhat  ahead  and 
behind  it  as  well.  The  key  idea  of  a  procedure  for  correcting  this,  pointed  out  by 
Nixon  (Refs.  14.15),  is  first  to  strain  the  coordinates  of  one  of  the  two  solutions  in 
such  a  fashion  that  the  shock  waves  align,  as  shown  In  the  upper  plot  on  the 
right  of  Figure  1 ,  and  then  determine  the  unit  perturbation.  Equivalently,  this 
can  be  considered  as  maintaining  the  shock  wave  location  invariant  during  the 
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perturbation  process,  and  assures  that  the  unit  perturbation  remains  small  both 
at  and  in  the  vicinity  of  the  shock  wave.  Obviously,  shock  points  are  only  one  of 
a  number  of  characteristic  high-gradient  locations  such  as  stagnation  points, 
maximum  suction  pressure  points,  etc.,  in  which  the  accuracy  of  the  perturbation 
soiution  can  degrade  rapidiy.  The  plots  in  the  lower  left  part  of  the  Figure  1 
indicate  such  a  situation  and  display  typical  supercritical  transonic  pressure 
distributions  which  contain  multiple  shocks  and  high-gradient  regions. 
Simultaneously  straining  at  all  these  locations  so  as  to  align  all  of  the 
characteristic  points  associated  with  each  shock  or  high-gradient  maxima  or 
minima,  as  indicated  in  the  lower  right  plot,  serves  to  minimize  the  unit 
perturbation  over  the  entire  domain  considered,  and  provides  the  key  to 
maximizing  the  range  of  validity  of  the  approximation  method.  Finally,  it  is 
important  to  recognize  that  while  the  approximation  method  replaces  the 
nonlinear  variation  between  two  nonlinear  base  solutions  with  a  linear  fit,  the 
resultant  approximation  solution  is  nonlinear  in  the  varied  parameter  because 
of  the  implicit  nature  of  coordinate  straining.  This  apparently  is  the  basis  of  the 
extended  range  of  accuracy  of  the  method. 

Because  the  method  of  strained  coordinates  is  known  to  be  nonapplicable  to 
certain  classes  of  perturbation  problems,  or  worse,  appear  to  be  applicable 
while  producing  incorrect  results  (Ref.  11).  the  point  arises  as  to  its  use  in  the 
present  context.  The  question  is  not  whether  the  method  will  produce  a 
uniformly  valid  solution  -  the  expansion  procedure  guarantees  that  -  but 
whether  the  solution  so  produced  may  be  incorrect  inherently.  Unfortunately, 
unlike  the  method  of  matched  asymptotic  expansions,  there  are  no  firm  rules 
which  guarantee  the  correctness  of  strained  coordinate  solutions.  There  are, 
nevertheless,  some  generally  reliable  guidelines.  The  method  of  strained 
coordinates  appears  to  always  succeed  when  the  singularity  predicted  by  the 
direct  unstrained  problem  actually  exists.  In  our  applications,  this  Is  always  the 
case,  since  we  identify  the  singularity  or  invariant  points  as  shock  points, 
stagnation  points,  and  other  physically  identifiable  points.  Furthermore,  we 
restrict  the  allowable  range  of  parametric  variation  such  that  the  neighboring 
calibration  and  predicted  flows  retain  these  same  points  and  create  no  new 
ones.  Therefore,  invariant  points  are  neither  lost  nor  generated  over  the 
solution  domain  of  interest.  These  considerations  effectively  Insure  that  the 
predicted  approximation  solutions  will  both  be  physically  correct  and, 
additionally,  will  not  violate  the  basic  straining  principle  (Ref.  1 1 )  regarding 
compounding  of  singularities. 


2.2  Previous  Applications  of  the  Approximation  Method 

At  this  point,  the  approximation  concept  based  on  the  ideas  discussed  above 
has  both  been  implemented  and  thoroughly  tested  In  a  wide  range  of  problems. 
In  Ref.  1,  several  candidate  ^proximation  methods  were  studied  and  the  most 
promising  method  was  identified.  Extensive  development  and  testing  of  that 
method  was  then  carried  out  in  Ref.  4  on  a  large  number  of  nonlinear  flow 
problems  involving  single-parameter  changes  of  a  variety  of  flow  and  geometric 
parameters.  Subcritical  and  supercritical  flows  past  isolated  airfoils  and 
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compressor  cascades  were  considered,  with  particular  emphasis  placed  on 
supercritical  transonic  flows  which  exhibited  large  surface  shock  movements 
over  the  parametric  range  studied.  Comparisons  of  the  approximation 
predictions  with  the  corresponding  "exact”  nonlinear  solutions  indicated  a 
remarkable  accuracy  and  range  of  validity  of  the  approximation  method.  For 
example.  Figure  2  from  Ref.  4  provides  a  comparison  of  results  illustrating  the 
remarkable  ability  of  the  approximation  method  to  predict  nonlinear  supercritical 
transonic  flows.  These  results  are  for  surface  pressures  obtained  from  full 
potential  solutions  and  represent  nonlifting  flows  past  several  NASA  four-digit, 
thickness-only  airfoils  at  Moo  »  0.820.  The  results  indicated  by  the  dotted  and 

dashed  lines  were  obtained  for  thickness  ratios  of  t  =  0.12  and  0.08, 
respectively.  Those  results  were  used  to  define  the  unit  perturbation  required 
by  the  approximation  method.  With  that  unit  perturbation  in  hand,  the 
approximation  method  was  then  employed  to  predict  surface  pressure  results 
for  thickness  ratios  x  =  (0.1 10.  0.105,  0.100,  0.095).  The  approximation  results, 
indicated  by  the  open  symbols,  v  are  then  compared  with  full  nonlinear  results 
obtained  by  running  the  full  potential  solver  at  those  thickness  ratios.  As  can  be 
seen,  the  results  are  essentially  identical,  in  particular,  in  the  region  of  the 
strong  shock. 

The  approximation  method  was  next  extended  in  Ref.  6  to  treat  simultaneous 
multiple-parameter  perturbations.  Extensive  testing  of  the  method 
demonstrated  remarkable  accuracy  and  range  of  validity  of  the  multiple- 
parameter  approximation  procedure  in  direct  correspondence  with  the  previous 
results  obtained  for  single-parameter  changes.  Additionally,  initial  applications 
of  the  multiple-parameter  approximation  method  combined  with  an  optimization 
procedure  were  also  made  to  several  two-dimensional  airfoil  design  problems. 
The  results  demonstrated  the  potential  of  the  approximation  method  for 
reducing  the  computational  work  in  certain  applications  by  an  order  of 
magnitude  with  no  degradation  in  accuracy.  Finally,  in  Ref.  7,  the 
approximation  method,  configured  in  a  form  suitable  for  predicting  an  arbitrary 
number  of  simultaneous  multiple  parameter  changes,  was  combined  with  the 
COPES/CONMIN  optimization  driver  (Ref.  16)  and  coupled  with  the  NASA 
TSONIC  full  potential  blade-to-blade  turbomachinery  solver  (Ref.  17).  A  series 
of  calculations  of  the  combined  code,  named  BLDOPT,  have  verified  the 
procedure,  demonstrated  the  accuracy  of  the  approximation-predicted  results, 
and  established  benchmark  guidelines  of  the  potential  for  computational 
savings  of  the  method  under  the  various  user  options  included  in  the  code.  In 
general,  the  approximation  method  was  found  to  be  capable  of  providing  an 
order  of  magnitude  reduction  in  computational  work  in  those  applications  which 
involved  essentially  subcritical  or  weakly  supercritical  turbomachinery  flows. 

2.3  Theoretical  Formulation:  Approximation  Method  Prediction  of  Surface 
Properties  on  3-D  Turbomachinery  Blades 

The  underlying  reason  for  the  remarkable  accuracy  of  the  approximation 
method  developed  in  this  study  lies  in  the  use  of  coordinate  straining  to  define 
the  unit  perturbation.  As  shown  in  Figure  1,  where  the  perturbation  between 
two  nonlinear  solution  states  is  displayed  graphically  as  the  shaded  area 
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between  the  base  and  the  strained  and  unstrained  calibration  solution, 
coordinate  straining  provides  the  ability  to  account  accurately  for  the 
displacement  of  a  multiple  number  of  discontinuities  and  maxima  of  hinh- 
gradient  regions  due  to  a  parameter  change.  This  enables  the  perturbation 
method  to  maintain  very  high  accuracy  in  regions  of  high  gradients  where  most 
perturbation  methods  commonly  fail,  and  to  maintain  that  accuracy  over  large 
parametric  ranges. 

In  what  follows,  we  provide  a  brief  account  of  the  theoretical  essentials  of  the 
strained-coordinate  perturbation  concept  as  configured  and  implemented  in  the 
present  design  application.  This  is  to  predict  simultaneous  multiple-parameter 
perturbation  flow  solutions  for  surface  properties  of  supercritical  wings  for  use  in 
optimized  wing  design.  The  flow  solutions  thus  considered  can  contain  a  total 
number  of  N  discontinuities  or  high-gradient  continuous  regions. 

To  proceed  with  the  theoretical  basis  of  the  approximation  method  as  applied  to 
simultaneous  multiple-parameter  perturbations  of  flows  containing  multiple 
shocks  or  high-gradient  regions,  consider  the  formulation  of  the  procedure  at 
the  full  potential  equation  level,  since  all  of  the  results  presented  here  are 
based  on  that  level.  Denote  the  operator  L  acting  on  the  full  velocity  potential  O 
as  that  which  results  in  the  three  dimensional  full-potential  equation  for  <!>,  i.e.. 


(1) 

If  we  now  expand  the  potential  in  terms  of  zero  and  higher-order  components  in 
order  to  account  for  the  variation  of  M  arbitrary  geometrical  or  flow  parameters  q) 
from  their  base  flow  values  qoj 


M 

0  =  0o+  X  qC>ij  +  ... 

j=1  * 

=  Poj  +  Aqj 


(2) 


and  then  insert  these  expansions  into  the  governing  Equation  (1),  expand  the 
result,  order  the  equations  into  zero  and  first-order  components,  and  make  the 

obvious  choice  of  expansion  parameters  ej  =  Aq)  we  obtain  the  following 
governing  equations  for  the  zero  and  M  first-order  components 


L[<Dol  =  0 

Lil<ti,l  +  ^M<l>ol  =  0 


(3) 


Here  Li  is  a  linear  operator  whose  coefficients  depend  on  zero-order  quantities 
and  3L[<I>oI/9qj  represents  a  "forcing"  term  due  to  the  qj*h  perturbation.  Actual 
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forms  of  Li  and  the  "forcing"  term  are  provided  in  Ref,  1  for  a  variety  of  flow  and 
geometry  parameter  perturbations  of  two-dimensional  turbomachines,  and  in 
Ref.  18  for  profile  shape  perturbations  of  an  isolated  airfoil.  An  important  point 
regarding  Equation  (3)  for  the  first-order  perturbations  is  that  these  equations 
represent  a  unit  perturbation  independent  of  the  actual  value  of  the  perturbation 

quantity  cj. 

Appropriate  account  of  the  movement  of  a  multiple  number  of  discontinuities 
and  maxima  of  high-gradient  regions  due  to  the  changes  in  the  parameters  qj  is 
now  accomplished  by  the  introduction  of  strained  coordinates  (s.t)  in  the  form 


M 

X  =  s  +  I  G  xi  (s.t) 

H 

M 

y=t+  I  ejyi(s.t) 
j=1 


(4) 


where 


N 

xi(s.t)  =  s+  X  6xi(t)xij(s.t) 
i=1 


N 

yi(s.t)=  I  5yjyij(s.t) 
i«i 


(5) 


and  EjSxi,  efiyj  represent  individual  x  and  y  displacements  due  to  perturbation  of 
the  qji^  parameter  of  the  N  strained  points,  and  xi|(s.t).  yi|(s.t)  are  straining 
functions  associated  with  each  of  the  N  strained  points.  For  the  applications 
considered  here,  we  have  assumed  that  all  discontinuities  such  as  shock  waves 
or  other  high-gradient  region  maxima  occur  essentially  normal  to  the  wing 
planform  so  that  only  the  (x.y)  coordinates  require  straining.  This  simplification 
is  not  strictly  necessary  and  could  be  relaxed  in  future  applications.  However, 
the  effect  of  this  assumption  on  the  prediction  of  surface  properties  via  the 
approximation  method  is  known  from  extensive  studies  of  the  two-dimensional 
case  to  be  of  higher  order  for  most  optimized  design  flow  situations  of 
aerodynamic  interest.  Introducing  the  strained  coordinate  Equations  (4)  and  (5) 
into  the  expansion  formulation  leaves  the  zero-order  result  in  Equation  (3) 
unchanged,  but  results  in  a  change  of  the  following  form  for  the  jth  perturbation 

Ll[4>l|l  +  L2,l<l)bl+ ^  L14tel  =  0  (6) 

Here  the  operators  are  understood  to  be  expressed  in  terms  of  the  strained  (s.t) 
coordinates,  and  the  additional  operator  L2j  arises  specifically  from 
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displacement  of  the  strained  points,  in  Refs.  15  and  18,  specific  expressions  for 
L2j  are  provided  for  selected  perturbations  involving  transonic  small- 
disturbance  and  full-potential  equation  formulations.  The  essential  point, 
however,  with  regard  to  perturbation  Equation  (6)  expressed  in  strained 
coordinates  is  that  it  remains  valid  as  before  for  a  unit  perturbation  and 
independent  of  ej. 

In  employing  the  approximation  method.  Equation  (6)  for  the  unit  perturbation 
is  solved  by  taking  the  difference  between  two  solutions  obtained  by  the  full 
nonlinear  procedure  after  appropriately  straining  the  coordinates.  If  we 
designate  the  solutions  for  some  arbitrary  dependent  flow  quantity  Q  as  base 
Qo  and  calibration  Qcj.  respectively,  of  the  varied  independent  parameter  qj,  we 
have  for  the  predicted  flow  at  some  new  parameter  value  qj  for  all  the  M 
parameters 


Q(x.y)  =  Qo(s.t)  + 


M 

I  qQij(s,t)  +  ... 


j=1 


(7) 


where 


'  q 

N 

Xj(s.t)  =  s  +  I  q  5xi(t)  xij  (s,t) 


N 

yj(t)  =  t  +  _I  q5y,yij(t) 


M 


x  =  s+  £ 
j=1 


S  (Xj(S.t)-S) 


M  £, 

y=t+  I  ^(yj(t)-t) 
j=i  ^ 


(8) 

(9) 

(10) 

(11) 

(12) 


q=qq-qoj  (13) 

q  =  qj-qoj  (14) 

We  note  that  in  order  to  determine  the  first-order  corrections  Qij(s,t),  we  require 
one  base  and  M  calibration  solutions  in  which  the  calibration  solutions  are 
determined  by  varying  each  of  the  M  arbitrary  independent  parameters  qj  by 
some  nominal  amount  from  the  base  flow  value  while  Keeping  the  other  fixed  at 
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their  base  values.  In  this  way,  the  first-order  corrections  Qii  can  be  determined 
from  Equation  (8)  where  Qq  is  defined  as  the  calibration  solution  corresponding 

to  changing  the  parameter  to  a  new  value  qcj.  Xj  is  the  strained  coordinate 

pertaining  to  the  Qcj  calibration  solution,  and  ej  >  qcj  -  qoj  represents  the  change 
in  the  qj  parameter  from  its  base  flow  value.  Thus, 


€j8xj(t)  =  (xf(t)-)^°(t))j 

(15) 

£j5Xi(t)=|  (xf(t)-Xj°(t))j 

(16) 

¥yi  =  (yVyTj)j 

(17) 

eiSyi=?  (yv^Ti)] 

(18) 

where  ej8xi(t)  given  in  Equation  (15)  represents  the  x  displacement  of  the  i*^ 
invariant  line  at  the  spanwise  t  location  in  the  p  calibration  solution  from  its 
base  flow  location  due  to  the  selected  change  ej  in  the  qj  parameter  given  by 
Equation  (13),  ej5xj(t)  given  in  Equation  (16)  represents  the  predicted  x 
displacement  of  the  i^^  invariant  line  at  the  spanwise  t  from  its  base  flow  location 
due  to  the  desired  change  Cj  in  the  qj  parameter  given  by  Equation  (1 4),  EjSyi 
given  in  Equation  (17)  represents  the  y  displacement  of  the  tip  of  the  i^^ 
invariant  line  in  the  calibration  solution  from  its  base  flow  location  due  to  the 
selected  change  ej  in  the  qj  parameter  given  by  Equation  (13).  and  ejSyi  given  in 
Equation  (18)  represents  the  predicted  y  displacement  of  the  tip  of  the  i*^ 
invariant  line  from  its  base  flow  location  due  to  the  desired  change  Ej  in  the  qj 
parameter  given  by  Equation  (14),  xij(s,t)  is  a  unit-order  straining  function 
having  the  property  that 


xii  (Xk(t).t)  «  I  J  (19) 

which  assures  alignment  of  the  1*^  invariant  line  at  the  t  spanwise  location 
between  the  base  and  calibration  solutions,  and  yi(t)  is  a  unit-order  straining 
function  having  the  property  that 


yii 


1 

0 


k  s  i 
k4i 


(20) 
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which  assures  alignment  of  the  tip  of  the  if**  invariant  line  between  the  base  and 
calibration  solutions. 

In  addition  to  the  single  conditions  given  by  Equation  (19)  and  (20)  on  the 
straining  functions,  it  may  be  convenient  or  necessary  to  impose  additional 
conditions  at  other  locations  along  the  contour.  For  example,  it  is  usually 
necessary  to  hold  invariant  the  end  points  along  the  contour,  as  well  as  to 
require  that  the  straining  vanish  in  a  particular  fashion  in  those  locations.  Ail  of 
these  conditions,  however,  do  not  serve  to  determine  the  straining  uniquely. 
The  nonuniqueness  of  the  straining,  nevertheless,  can  often  be  turned  to 
advantage,  either  by  selecting  particularly  simple  classes  of  straining  functions 
or  by  requiring  the  straining  to  satisfy  further  constraints  convenient  for  a 
particular  application. 

The  fact  of  nonuniqueness  of  straining  function,  however,  raises  a  further 
question  of  the  dependence  of  the  final  approximation-predicted  result  on 
choice  of  straining  function.  An  initial  example  of  the  effect  of  employing  two 
different  straining  functions  for  a  strongly  supercritical  two-dimensional  flow  was 
provided  in  Ref.  15,  and  in  Ref.  4  a  detailed  examination  was  made  of  the 
dependence  of  approximation  results  on  several  classes  of  different  straining 
functions.  Although  it  can  be  demonstrated  (Ref.  15)  that  the  final 
approximation-predicted  result  obtained  when  employing  strained  coordinates 
is  formally  independent  of  the  particular  straining  function  used  -  provided  that 
the  straining  function  moves  the  invariant  points  to  the  proper  locations  -  the 
results  of  Ref.  4  demonstrate  that,  under  certain  conditions,  particular  classes  of 
straining  functions  can  induce  spurious  approximation  results.  The  underlying 
reason  is  that,  while  the  approximation-predicted  results  at  and  in  the  vicinity  of 
invariant  points  are  independent  of  the  choice  of  straining  function  (provided 
invariant  point  locations  are  preserved),  some  classes  of  straining  functions 
have  the  undesirable  property  of  producing  unwanted  straining  in  certain 
regions  removed  from  the  invariant  points.  The  correction  for  this  deficiency, 
which  was  found  in  Ref.  4  and  has  proven  effective  in  all  case  studies 
undertaken,  is  to  employ  linear  piecewise-continuous  straining  functions.  This 
both  preserves  the  accuracy  of  the  approximation  results  in  the  vicinity  of  the 
invariant  points,  and  introduces  no  excessive  straining  in  regions  removed  from 
those  locations. 


For  linear  piecewise-continuous  straining  functions,  the  functional  forms  of  the 
straining  can  be  compactly  written.  For  the  x  displacement  we  have 


Xj  (s,t)  =  s  +  ^ 


xw(t) 


'41 


-  X; 


(xf  (t)  •  X°(t))j 


s-)^?(t) 

X4i(t)-x°(t) 


(x|;i(t)-Xi°i(t))j 


•H(Xj°^(t)-s)*H(s-Xj°(t)) 


(21) 
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where  H  denotes  the  Heaviside  step  function.  As  discussed  above,  it  is  usually 
necessary  to  hold  invariant  both  of  the  end  points  along  the  contour  in  addition 
to  the  points  corresponding  to  discontinuities  or  high-gradient  maxima. 
Consequently,  for  the  results  reported  here,  the  array  of  x  invariant  points  in  the 
base  and  calibration  solutions  have  been  taken  as 

x°(t)  =  {0,  x°(t),  4(0 .  4(0- 


xfj(t)  =  {0,  x®|(t),  x|(t), ....  x®j(t).  1} 


(22) 


where  the  contour  length  at  the  spanwise  location  t  has  been  normalized  to 
unity  and  where  n  is  the  number  of  invariant  points  along  the  contour  exclusive 
of  the  end  points. 

Similarly,  for  the  y  displacements  of  the  tips  of  the  continuity  lines  we  have 


yj(t)  =  t  +  ^ 


(yv»T|)i*-5 


1+1 


o 


,  c 

•(VTi 


1+1  yVi^i 


•  H(y?.^^  - 1)  •  H(t  -  y?.) 


(23) 


where  the  spanwise  locations  of  the  tips  of  the  discontinuity  lines  in  the  base 
and  calibration  solutions  have  been  taken  as 


0  ,0  0  0  0  . 
yjj  =  (vti'  Via*  yja*  ’  yin^ 

yiij  =  (vti’  ^12’  yja . Vin^j 


(24) 
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3.  RESULTS 


3.1  Development  of  Multiple  invariant  Point  Characterization  Procedure  for 
Surface  Shock  Waves 

Development  of  an  accurate  procedure  to  enable  the  definition  and  subsequent 
tracking  of  the  invariant  points  associated  with  high-gradient  regions,  in 
particular,  shock  waves  is  key  to  the  accuracy  of  the  ^proximation  method. 
Both  current  and  previous  testing  of  the  approximation  method  has 
demonstrated  the  advantages  of  employing  a  multiple  invariant  point 
representation  of  shock  waves.  Results  for  the  approximation  prediction  of 
highly  nonlinear  solutions  via  the  coordinate  straining  procedure  employed  in 
the  method  clearly  show  that  the  accuracy  of  the  predicted  results  are  most 
sensitive  in  the  vicinity  of  shock  waves  or  other  high  gradient  compressive 
regions.  This  occurs  since  in  those  regions  any  misalignment  of  the  base  and 
calibration  solutions  employed  to  determine  the  unit  perturbations  from  which 
the  approximate  solutions  are  obtained  becomes  highly  magnified  in  the 
approximation  predicted  results  due  to  the  presence  of  the  large  gradients  in 
these  regions.  Hence,  a  single  invariant  point  representation  of  a  shock  is  quite 
accurate  for  normal  or  near-normal  shock  waves.  However,  for  shocks  that  are 
characterized  by  more  gradual  gradients  so  that  the  flow  variation  is  that  of  a 
highly  sloped  but  not  vertical  line,  predictions  based  on  the  single  invariant 
point  representation  can  degrade  rapidly  for  even  moderate  parameter 
extrapolation  or  interpolations  from  the  base  and  calibration  values.  A  method 
for  alleviating  this  degradation  in  accuracy  is  to  introduce  additional  invariant 
points  to  characterize  the  shock  wave.  This  device  acts  to  improve  the 
predicted  approximate  solution  in  those  high  gradient  regions  by  requiring  the 
approximate  method  to  track  more  than  one  topological  feature  of  the  shock. 
Shock  waves  that  are  typically  generated  by  most  of  the  currently  existing 
Navier-Stokes,  Euler,  and  full  potential  flow  algorithms  are  usually 
characterized  by  some  precursor  compression  region  upstream  of  the  shock,  a 
high-gradient  central  shock  region,  and  a  post  expansion  region,  such  as 
illustrated  in  the  sketch  below. 


Illustration  of  Multiple  Invariant  Point  Characterization 
for  Surface  Shock  Wave  Topology  Tracking 
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For  our  purposes,  a  sufficiently  general  way  to  characterize  such  a  shock  would 
be  to  use  four  Invariant  points  that  would  correspond  to:  (1)  pre-shock  minimum 
pressure  point,  (2)  maximum  gradient  or  Cp*  point,  (3)  post-shock  maximum 

pressure  point,  and  (4)  post-shock  minimum  expansion  pressure  point.  While 
such  a  characterization  of  a  shock  is  quite  general,  it  is  not  necessarily  the  most 
advantageous  choice  to  employ  in  all  cases.  The  reason  is  that  the  local 
topology  of  shock  waves  in  many  transonic  flows,  even  over  relatively  small 
parameter  ranges,  can  change  so  as  to  make  it  difficult  to  continuously  identify 
all  the  invariant  points  described  above.  For  example,  in  many  transonic  flows  it 
is  common  for  a  change  in  the  post-shock  expansion  topology  to  occur  during  a 
parameter  variation  whereby  that  region  loses  its  post  expansion  character. 
The  remedy  for  this  is  to  use  only  the  minimum  number  of  invariant  points 
needed  to  track  a  topological  phenomenon  over  a  parameter  range  of  interest. 

With  regard  to  shock  waves,  we  have  concluded  that  for  the  class  of  flows  of 
primary  interest  to  this  study  that  it  is  safe  to  ignore  the  post  shock  expansion 
region  (i.e.  point  4)  and  confine  the  multiple  invariant  point  characterization  to 
the  compressive  gradient  portion  of  the  shock,  i.e.  between  points  (1 )  and  (3) 
illustrated  above.  We  have  done  this,  and  have  found  from  numerical 
experiments  that  the  post  expansion  region  remains  satisfactorily  predicted. 
Additionally,  we  have  found  that  employing  a  third  invariant  point  located  at  the 
maximum  gradient  or  Cp*  location  in  the  middle  of  the  shock  (point  2)  adds  no 
significant  improvement  in  accuracy.  Hence,  for  the  transonic  flows  of  concern 
to  this  study,  a  two  invariant  point  characterization  of  shock  waves  is  sufficient 
and  was  adopted  as  the  standard  model. 

In  order  to  track  the  two  invariant  shock  points,  it  is  necessary  to  establish  an 
appropriate  criteria  to  continuously  identify  them.  Ideally,  the  criteria  should  be 
such  so  as  to  allow  unambiguous  identification  of  the  invariant  points 
throughout  the  anticipated  range  of  parameter  variation.  In  order  to  establish 
these  criteria  and  because  these  shock  wave  invariant  points  play  such  a 
critical  role  in  controlling  the  overall  accuracy  of  the  approximation  method,  it 
was  necessary  to  carry  out  extensive  numerical  experimentation  to  identify  the 
best  criteria  to  set  these  points.  As  a  result  of  this  test,  we  have  found  the 
following  strategy  to  identify  and  track  the  two  shock  wave  invariant  points  the 
best  of  a  variety  of  schemes. 

Selection  of  the  invariant  point  at  the  beginning  of  the  shock  was  found  to  be 
best  set  by  choosing  the  point  where  dCp/dx  =  0.  We  investigated  other 
alternatives,  such  as  placing  the  invariant  point  where  the  shock  slope  attains  a 
certain  value  (eg.  dCp/dx  s  -i)  or  placement  of  the  invariant  point  at  a  certain 
location  upstream  of  the  point  where  the  pressure  distribution  goes  through 
sonic  (Cp  »  Cp* ).  The  result  of  these  studies  has  shown  that  the  choice  of  the 
point  where  dCp/dx  s  0  is  the  best. 

Selection  of  the  invariant  point  at  the  foot  or  end  of  the  shock  is  more  direct  in 
the  sense  that  the  maximum  pressure  point  behind  the  shock  is  a  natural  and 
obvious  choice  for  termination  of  the  shock  region.  However,  while  the  choice 
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of  the  point  is  clear,  we  have  found  from  the  results  of  our  testing  that  in  many 
instances  the  approximation  predictions  in  the  vicinity  of  the  shock  foot  degrade 
rapidly.  Furthermore,  we  found  that  significantly  improved  results  in  that  region 
could  be  achieved  by  modification  of  the  location  of  the  shock  foot  point.  These 
cases  occur  when,  due  to  local  numerical  grid  spacing  of  the  flow  solver  or 
other  factors,  the  maximum  pressure  point  behind  a  shock  in  the  base  and 
calibration  solutions  happens  to  lie  just  off  the  maximum  gradient  line  of  the 
shock,  as  illustrated  in  the  sketch  below  on  the  right,  rather  than  on  the 
maximum  gradient  line  as  normally  occurs  as  shown  in  the  figure  on  the  left. 


Shock  Foot  Point  On 
Max  Gradient  Line 


Shock  Foot  Point  Off 
Max  Gradient  Line 


The  remedy  for  this  is  to  locally  sharpen  the  shock  prior  to  determining  the  unit 
perturbation  distributions  by  relocating  the  shock  foot  on  the  maximum  gradient 
line.  We  have  implemented  this  into  the  approximation  method,  and  have  found 
that  this  results  in  significantly  improved  results  in  the  predicted  approximate 
solutions  for  both  the  location  of  the  shock  foot  as  well  as  the  pressure 
distribution  in  the  near  vicinity  of  the  shock  foot. 

Finally,  even  with  all  the  careful  analysis  involved  in  the  selection  of  the  two 
invariant  points  characterizing  the  shock,  it  is  inevitable  in  many  cases  involving 
multiple  parameter  perturbations  that  there  will  often  occur  some  slight 
misalignment  in  the  shock  regions  between  the  base  and  the  various  calibration 
solutions.  This  results  in  inaccuracies  in  the  unit  perturbations  and.  as  a 
consequence,  in  the  predicted  approximation  solutions.  Particularly  in 
situations  involving  extreme  parameter  extrapolation  or  interpolation,  the 
predicted  shock  region  often  contains  a  point  or  points  for  which  the  results  are 
clearly  spurious.  While  the  most  direct  means  of  eliminating  this  would  be  to 
employ  finer  grids  for  the  computation  of  the  base  and  calibration  flow  solutions, 
this  is  often  not  desirable  or  possible  because  of  the  increased  computer 
storage  and/or  CPU  costs.  Consequently,  what  is  required  is  a  procedure  for 
smoothing  the  predicted  values  at  those  points  while  simultaneously  retaining 
the  overall  integrity  of  the  approximation  solution.  We  have  implementated  a 
combination  of  several  schemes  that  ensure  the  appropriate  behavior  of  the 
predicted  solution  in  the  shock  region.  To  accomplish  this,  it  was  necessary  to 
employ  more  than  one  smoothing  procedure  in  order  to  allow  for  the  correction 
of  from  one  to  a  number  of  points  in  the  shock  region.  For  example,  we  have 
found  it  best  to  employ  a  linear  smoothing  correction  in  order  to  correct  for  one 
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or  two  points,  while  for  three  or  more  points  it  was  found  to  be  necessary  to 
employ  either  quadratic  or  cubic  least  squares  smoothing.  The  smoothing 
procedures  developed  ensure  that  monotonidty  is  maintained  in  both  ordinate 
and  slope  throughout  the  predicted  region. 

With  these  procedures  incorporated  into  the  approximation  method,  we  have 
performed  extensive  testing  to  verify  them  on  difficult  one,  two.  and  three 
parameter  perturbations  involving  highly  nonlinear  supercritical  transonic  flows. 
Figures  3  through  6  provide  some  examples  of  the  remarkable  predictions  that 
the  approximation  method  is  capable  of. 

Figures  3a.-g.  presents  the  results  of  the  ability  of  the  approximation  method  to 
predict  extreme  solution  extrapolation  results  when  using  two  very  closely 
spaced  base  and  calibration  solutions.  The  base  and  calibration  solutions 
shown  and  used  to  determine  the  unit  perturbations  needed  for  the 
approximation  results  are  for  highly  supercntical  flows  past  a  symmetric  NACA 
OOXX  series  blade  profile  at  an  oncoming  Mach  number  of  0.82  and  zero  angle 
of  attack  with  thickness  ratios  t/c  =  0.115  and  0.120  for  the  base  and  calibration 
solutions,  respectively.  The  predicted  results  indicated  by  open  symbols  are  for 
new  thickness  ratios  t/c  =  {0.110,  0.105,  0.100,  0.095,  0.090,  0.085,  0.080}  and 
are  shown  in  Figures  3a.  to  3g.  Those  results  are  meant  to  be  compared  with 
the  solid  lines  which  represent  the  exact  solutions  obtained  from  running  the 
flow  solver  at  the  new  thickness  ratios.  We  believe  these  results  are 
remarkable  in  that  even  though  the  base  and  calibration  solutions  were 
purposely  li^orly  selected  for  spanning  the  range  of  thickness  ratios  shown,  the 
approximation  method  is  able  to  do  an  excellent  job  of  tracking  the  shock  for 
extrapolations  beyond  400%  of  the  parameter  difference  between  the  base  and 
calibration  solutions.  This  is  far  beyond  what  any  linear  approximation  theory  is 
capable  of  achieving. 

The  corresponding  results  shown  in  Figures  4a.-g.  are  for  a  more  reasonable 
choice  of  base  and  calibration  solutions.  For  these  results,  the  base  and 
calibration  thickness  ratios  were  t/c  s{0.900.0.110}.  respectively.  The  predicted 
results  are  for  t/c  =  (0.120,  0.115,  0.105,  0.100,  0.095,  0.085,  0.080},  and 
represent  both  solution  interpolation  and  extrapolation.  The  approximation 
predictions  agree  extremely  well  with  the  exact  results  over  the  entire  parameter 
range  displayed. 

In  Figures  5a.-c.  we  provide  results  demonstrating  the  capability  of  the 
approximation  method  to  simultaneously  predict  both  upper  and  lower  surface 
pressures.  The  results  are  for  an  angle  of  attack  perturbation  of  highly 
supercritical  flows  past  a  NACA  0012  blade  profile  with  oncoming  Mach  number 
M  s  0.75.  The  base  and  calibration  angles  of  attack  were  {1.00,  2.00}  degrees, 
respectively,  with  the  predicted  angles  at  {2.50, 1.50,  0.50}.  The  approximation 
results  indicated  are  again  in  good  agreement  with  the  exact  solutions  over  the 
entire  parameter  range  and  on  both  surfaces  of  the  blade. 

Finally,  in  Figures  6a.-e.  we  provide  the  results  of  a  simultaneous  two- 
parameter  perturbation  of  strongly  nonlinear  flows.  The  results  shown  are  for 
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an  oncoming  Mach  number  and  thickness  ratio  perturbation  of  a  symmetric 
NACA  OOXX  blade  profile  at  zero  angle  of  attack.  The  base  and  calibration 
Mach  number  and  thickness  ratio  are  1^  »  0.62,  (t/c)b  «  0.120,  Me  «  0.80,  (t/c)c 
=  0.110.  The  predicted  results  are  for  {M,  t/c}  >  {0.84,  0.115),  {0.83,  0.115), 
{0.81,  0.115),  {0.79,  0.115),  {0.78,  0.115).  The  agreement  between  the 
approximation  predicted  and  exact  nonlinear  results  are  again  excellent.  We 
have  performed  many  additional  comparative  caculations  with  the 
approximation  method  and  have  obtained  similar  good  results  as  those  shown 
above. 


3.2  Multiple  Invariant  Point  Characterization  and  Invariant  Point 
Preprocessing  Procedure 

The  first  step  in  developing  the  approximation  method  solution  is  the 
determination  of  the  particular  invariant  points  to  be  employed  to  characterize 
the  flows  under  consideration.  For  the  typical  nonlinear  flows  pertinent  to  this 
study,  these  invariant  points  can  be  characterized  as  follows:  maximum 
pressure  or  stagnation  point,  minimum  suction  pressure  and  associated 
recovery  points  near  the  nose  on  both  the  upper  and  lower  surfaces.  Cp*  points 

denoting  shock  waves  and  their  associated  local  minimum  and  maximum 
pressure  points  that  characterize  the  initiation  and  termination  of  the  shock 
regions.  Figure  7  illustrates  a  3-D  pressure  distribution  characteristic  of  the 
ones  of  interest  to  the  present  study.  Shown  in  that  figure  are  the  chordwise 
surface  pressure  distributions  on  both  the  upper  and  lower  surfaces  of  a  3-D 
blade  at  17  spanwise  stations  across  the  blade  from  root  to  tip.  In  locating  the 
invariant  points  for  this  example,  we  begin  at  the  trailing  edge  and  consider  the 
trailing  edge  points  for  both  the  upper  and  lower  surface  fixed.  Thus,  we 
proceed  from  the  trailing  edge  point  on  the  lower  surface  and  move  upstream 
searching  for  sonic  Cp*  points  and  the  stagnation  point.  If  a  Cp*  point  is  found 
on  the  lower  surface,  we  then  proceed  to  locate  a  local  dCp/dx  »  0  point 
upstream  of  the  shock  and,  if  appropriate,  search  for  a  local  maximum  pressure 
point  just  downstream  of  the  Cp*  point.  After  locating  all  the  Cp*  points  and 
their  associated  maximum  and  minimum  points  on  the  lower  surface,  we  locate 
the  maximum  pressure  or  stagnation  point  near  the  nose.  We  then  locate  the 
minimum  suction  pressure  point  on  the  lower  surface  Just  downstream  of  the 
stagnation  point,  and  the  dCp/dx  s  0  point  just  downstream  of  the  tower  surface 
minimum  suction  pressure  point.  Next,  the  minimum  suction  pressure  point  on 
the  upper  surface  located  just  downstream  of  the  stagnation  point  is  found,  as 
well  as  the  associated  dC^^dx  >  0  points  just  upstream  of  those  points  and  the 

maximum  post-shock  pressure  points  just  downstream  of  those  Cp*  points.  As 

an  example,  for  the  sample  case  shown  in  Figure  7,  nine  invariant  points  were 
located  according  to  the  above  criteria.  These  invariant  points  on  each  of  the 
17  chordwise  pressure  distributions  are  indicated  by  a  (-i-)  sign.  They 
correspond  to  a  Cp*  and  assodated  dCp/dx  s  0  point  on  the  tower  surface  near 
the  nose,  the  stagnation  point,  the  minimum  suction  pressure  and  associated 
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dCp/dx  >  0  point  on  the  upper  surface  near  the  nose,  and  the  maximum  post¬ 
shock  and  associated  upstream  dCp/dx  >  0  point  for  the  upper  surface  shock. 

The  next  step  in  proceeding  with  the  approximation  method  solution  is  to 
preprocess  these  numerically  computed  pressure  distributions  at  each 
spanwise  location  so  as  to  smooth  them  in  the  vicinity  of  the  shock  waves  and 
minimum  suction  pressure  locations.  The  result  of  this  preprocessing  for  the 
example  shown  in  Figure  7  is  illustrated  in  Figure  8.  Comparison  of  the  results 
in  Figure  8  with  those  of  Rgure  7  shows  that  as  a  result  of  the  preprocessing  a 
distinct  sharpening  of  the  shock  in  the  vicinity  of  its  post-shock  maximum 
pressure  has  been  achieved. 

The  final  step  of  the  preprocessing  procedure  is  to  relocate  all  these  invariant 
points  in  each  of  the  chordwise  pressure  distributions  so  as  to  ensure  a  globally 
continuous  variation  in  the  spanwise  coordinate  direction  of  the  various 
selected  points.  The  result  of  this  final  step  of  the  smoothing  procedure  is 
illustrated  in  Figure  9.  Comparison  of  the  relocated  invariant  points  with  the 
original  invariant  points  shown  in  Figure  7  clearly  displays  a  smooth  global 
variation  is  now  present  in  the  locations  of  both  the  dCp/dx  »  0  points 
associated  with  the  upper  surface  minimum  suction  pressure  point,  and  also 
with  the  upper  surface  post-shock  maximum  pressure  points. 

In  order  to  demonstrate  how  accurately  the  approximation  method  enhanced  by 
the  new  preprocessing  scheme  can  predict  strongly  nonlinear  3-D  flows,  we 
have  selected  two  different  3-D  test  problems.  The  first  involves  an  overall 
thickness  ratio  change  of  an  isolated  3-D  blade  having  an  ONERA  M6  profile. 
To  establish  the  parameter  range  of  thickness  ratios  over  which  the 
approximation  method  would  be  tested,  we  considered  changes  in  the  overall 
thickness  ratio  of  the  blade  from  80%  to  120%  in  steps  of  5%  from  its  nominal 
value,  i.e.  (Vc)  =  {0.80,  0.85,  0.90,  0.95, 1.00, 1.05, 1.10, 1.15, 1.20}  so  that  nine 
3-D  solutions  were  determined.  These  were  computed  at  an  oncoming  Mach 
number  of  M  s  0.86  and  angle  of  attack  of  3.06  degrees.  After  application  of  the 
preprocessing  precedures  described  above,  the  solution  for  the  nominal  value 
of  thickness  ratio  (t/c)ff,ax  «  1.00  is  illustrated  in  Figure  10  where  the  chordwise 
upper  and  lower  surface  pressure  distributions  at  each  of  the  17  spanwise 
locations  of  the  3-D  grid  used  in  the  solution  determination  are  displayed. 
Figures  11a.-c.  provide  comparisons  of  the  approximation  method  at  the  root 
chord  station  (y/s  =  0.0)  when  the  base  and  calibration  solutions  are  taken  at 
(l/c)max  {0.90,  1.10),  respectively,  and  the  predicted  results  are  given  at 

(^/c)max  (0.80,  1.00,  1.20).  The  comparisons  between  the  approximation 
meth(^  and  the  exact  results  are  outstanding.  Rgures  12a.-c.  provide  similar 
results  at  the  outermost  spanwise  station  (y/s  »  0.969).  Again  the  comparisons 
are  remarkable. 

The  second  3-D  example  problem  involves  a  significantly  more  difficult  test  of 
the  approximation  method  due  to  the  fact  that  we  selected  it  to  involve  an 
unusually  large  topology  change  in  the  3-D  solutions  over  the  parameter  range 
selected.  This  problem  involves  an  oncoming  Mach  number  parameter  change 
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over  the  same  3-D  blade  geometry  studied  in  the  first  example,  but  with  the 
Mach  number  changing  from  M  »  0.880  to  0.720  in  steps  of  0.02.  Three  of  the 
nine  3-D  solutions  at  the  beginning,  middle,  and  end  of  the  parameter  range 
after  invariant  point  preprocessing  are  displayed  in  Figures  13a.-c.  The  large 
topology  change  is  evident  as  we  go  from  the  strongly  supercritical  pressure 
distribution  at  M  »  0.880  shown  in  Figure  13a..  which  has  an  extremely  strong 
upper  surface  shock,  to  the  subcritical  distribution  at  M  s  0.720  shown  in  Figure 
13c..  which  has  no  upper  surface  shock  but  has  an  extremely  pronounced  and 
compressed  upper  surface  minimum  suction  pressure  region  near  the  nose. 
Although  such  a  major  topology  change  is  normally  considered  far  beyond  the 
assumptions  under  which  the  approximation  method  is  valid,  we  have  begun 
development  of  an  extension  of  the  method  to  allow  it  to  at  least  approximate 
solutions  through  these  topology  changes  without  further  recourse  to  yet 
another  new  set  of  3-D  base  and  calibration  solutions  that  would  cover  the  new 
topology  change.  This  procedure  involves  development  of  a  criteria  to  handle 
invariant  point  coalescence  whereby  two  invariant  points  move  toward  one 
another  and  eventually  meet.  When  that  occurs,  a  new  criteria  is  necessary  to 
implement  the  approximation  method  in  order  to  prevent  the  two  invariant  points 
from  crossing  and  producing  topologically  spurious  results.  In  the  planned 
optimization  applications  planned  for  the  approximation  method,  we  believe  it  to 
be  extremely  advantageous  if  the  method  can  continue  to  perform  reasonably 
well  and  in  an  automatic  fashion  even  under  such  extended  circumstances. 
Some  results  of  the  approximation  method  predictions  under  such  conditions 
are  provided  in  Rgures  14a.-d.  These  results  are  at  the  most  extreme  spanwise 
location  (y/s  «  0.969)  along  the  blade,  and  are  for  a  choice  of  base  and 
calibration  Mach  numbers  of  M  ==  {0.82.  0.86}.  with  the  predictions  given  at  M  > 
(0.88.  0.84,  0.80.  0.72}.  We  clearly  see  the  disparate  differences  between  the 
base  and  calibration  solutions  in  these  plots.  Nevertheless  with  regard  to  the 
comparisons,  the  approximation  method  is  able  to  do  a  very  reasonable  job  of 
predicting  the  pressure  distribution  change  from  strongly  supercritical  as  shown 
in  Figure  14a.  to  quite  subcritical  as  shown  in  Figure  14d.  An  additional 
procedure  to  eliminate  the  spurious  points  evident  behind  the  upper  surface 
minimum  suction  pressure  peak,  and  also  to  improve  the  approximation 
prediction  behind  the  maximum  suction  pressure  point  is  needed. 


3.3  Coupling  of  3-D  Approximation  Method  and  3-D  Flow  Solver  With  Design 
Optimization  Procedures 

3.3.1  Selection  of  Design  Optimization  Procedure 

A  review  was  conducted  of  all  methods  and  applications  thereof  that  have  been 
developed  and  applied  to  the  general  class  of  aerodynamic  design  optimization 
problems  that  are  of  relevance  to  this  study.  The  concept  of  optimizing  about  an 
aerodynamic  rather  than  structurally  related  design  objective  is  relatively  new 
(Refs.  19,20).  Consequently,  with  one  exception,  comparative  studies  of  the 
performance  and  reliability  of  different  optimization  algorithms  for  this  class  of 
optimization  problems  have  not  been  earned  out.  The  result  has  been  that 
essentially  all  of  the  current  applications  involving  aerodynamic  optimization 
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studies  (eg.  Refs.  21-25)  have  involved  use  of  the  CONMIN  conjugate  gradient 
optimization  method  (Ref.  26),  a  procedure  originally  designed  for  the 
conceptually  different  problem  of  optimized  structural  design.  The  CONMIN 
method,  as  well  as  the  more  generalized  COPES  procedure  (Ref.  16)  which 
includes  CONMIN  as  the  basic  optimization  module  but  contains  additional 
modules  that  provide  some  convenient  features  for  carrying  out  certain 
constrained  minimization  problems  involving  side  constraints  and  design 
variable  bounding,  have  now  been  applied  to  a  variety  of  sensitive 
aerodynamic  design  optimization  problems.  These  have  involved  the  optimized 
design  of  airfoils  (Refs.  19,20),  turbomachinery  blades  (Ref.  7)  and  wings 
(Refs.  23-25)  at  highly  nonlinear  supercritical  transonic  flow  conditions.  For  the 
structural  optimization  problems  for  which  the  CONMIN  procedure  was 
originally  developed,  the  method  has  generally  performed  quite  satisfactorily. 
However,  for  the  new  class  of  aerodynamic  optimization  problems,  the  CONMIN 
procedure  has  demonstrated  a  number  of  deficiencies.  This  is  not  altogether 
surprising  since  structural  design  optimization  problems  typically  involve  large 
numbers  (10  -  1000)  of  design  variables  and  side  constraints  but  rapid 
computation  of  the  governing  internal  structural  element  solutions.  In  contrast, 
aerodynamic  optimization  problems  typically  involve  relatively  few  design 
variables  and  side  constraints  but  slow  and  expensive  computations  of 
relatively  sensitive  governing  flow  solutions. 

We  have  reviewed  all  the  published  applications  involving  aerodynamic 
optimization  problems  of  the  general  class  of  interest  here  in  which  the 
CONMIN  optimization  procedure  has  been  employed.  These  investigations 
include  our  own  substantial  experiences  with  the  CONMIN  and 
COPES/CONMIN  procedures  in  which  we  designed  a  pilot  code  for  optimizing 
the  design  of  biade-to-blade  profiles  of  turbomschinery  blades  (Ref.  7).  In 
assessing  the  performance  of  the  CONMIN  optimizer,  which  employs  the 
Fletcher-Reeves  conjugate  gradient  algorithm  as  the  driver,  all  of  the  reported 
investigations  appear  to  have  encountered  similar  deficiencies  with  the  method. 
The  most  significant  of  those  is  the  tendency  of  the  conjugate  gradient  algorithm 
to  become  focused  or  stuck  at  a  local  minimum,  and  not  be  able  to  discern  that 
this  has  occurred  and  proceed  from  there  toward  a  more  global  minimum.  In 
the  worst  cases,  the  method  does  not  move  to  any  significant  degree  away  from 
the  original  design.  The  basis  of  this  problem  lies  with  the  manner  in  which  the 
conjugate  gradient  method  performs  its  optimization  search.  At  each  iteration  of 
the  techriique,  the  gradient  of  the  objective  function  with  respect  to  each  design 
variable  is  calculated  and  a  line  search  is  performed  in  the  direction  given  by  a 
certain  linear  combination  of  the  current  gradient  and  the  last  search  direction. 
The  objective  function  gradients  are  always  calculated  by  one-sided  forward 
differences  with  no  provision  for  central  differences.  The  step  size  employed 
can  be  scaled  by  the  size  of  the  associated  variables,  but  is  independent  of  the 
value  of  the  objective  function,  its  precision  or  its  derivatives.  What  this 
ultimately  mearis  is  that  any  variable  which  increases  the  objective  function 
when  stepped  in  one  direction  is  assumed  to  reduce  the  objective  function 
when  moved  in  the  opposite  direction.  This  is  not  true  when  the  objective  is 
near  a  local  minimum  of  one  variable.  It  Is  in  these  regions  where  the  method 
can  fail  disastrously.  What  can  occur  near  these  local  minima  is  that  as  the 
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gradient  alternates  in  sign  about  that  point,  the  stepsize  is  reduced  by  the 
optimizer  in  hopes  of  resolving  the  problem.  With  the  reduced  stepsize,  the 
solution  process  becomes  unable  to  move  away  from  that  location.  This  result 
is  exactly  what  has  been  observed  with  the  CONMIN  conjugate  gradient 
method. 

One  way  to  try  to  alleviate  this  problem  would  be  to  employ  central  differencing 
in  the  gradient  determinations.  While  this  would  be  substantially  more 
expensive  than  employing  one-sided  differences,  it  would  not  have  to  be 
implemented  at  all  times,  but  only  when  the  optimization  search  is  in  the  vicinity 
of  a  local  minimum.  We  have  investigated  implementing  this  possibility  by 
examining  the  combined  CONMIN/TSONIC  blade  optimization  code  that  we 
previously  developed  (Ref.  7)  to  determine  whether  central  differencing  could 
be  easily  introduced  into  the  CONMIN  optimization  procedure.  It  appears  that 
implementing  a  central  difference  procedure  into  the  method  cannot  be  easily 
done,  but  would  require  a  significant  rewrite  of  the  optimization  code.  This 
would  be  an  effort  beyond  the  scope  of  the  present  investigation. 

The  one  tested  alternative  to  the  COPES/CONMIN  conjugate  gradient 
procedure  for  aerodynamic  optimization  applications  is  the  quasi-Newtonian 
procedure.  That  procedure  has  now  been  implemented  into  an  optimization 
driver  called  QNMDIF  (Ref.  9),  and  configured  so  as  to  avoid  several  of  the  more 
serious  deficiencies  inherent  in  the  CONMIN  driver.  The  QNMDIF  driver  is 
based  on  the  quasi-Newtonian  minimization  algorithm  as  orginally  developed 
by  Gill,  et.  al.  (Refs.  27,  28).  A  number  of  improvements  have  been 
implemented  to  maintain  stable  convergence  in  the  presence  of  roundoff  errors, 
ill-conditioning,  or  occasional  small  discontinuities  in  the  objective  function.  For 
example,  the  method  will  switch  from  fonward  to  central  differences  if  the  search 
procedure  fails  to  produce  a  significantly  better  point,  and  then  back  to  forward 
differences  if  an  improved  rate  of  progress  occurs  later.  A  comparative  series  of 
systematic  test  cases  employing  both  the  CONMIN  and  QNMDIF  solvers  have 
been  carried  out  for  a  number  of  difficult  optimization  problems  (Ref.  9).  The 
results  have  demonstrated  the  significantly  better  performance  of  the  quasi- 
Newton  algorithm.  The  QNMDIF  procedure,  originally  developed  for  2-D 
applications,  has  now  been  combined  with  the  3-D  TWING  full  potential  flow 
solver  into  &.n  efficient  optimization  so'ver  for  supercritical  3-D  wing  design  (Ref. 
8). 

On  the  basis  of  the  demonstrated  superiority  of  this  method  over  the  CONMIN 
conjugate  gradient  procedure,  we  have  opted  to  employ  the  QNMDIF  quasi- 
Newton  procedure  as  the  optimization  driver  of  choice  for  the  present 
investigation.  The  integration  of  the  3-D  approximation  method,  as  currently 
developed,  with  the  combined  3-D  flow  field  solver  (TWING)  and  quasi-Newton 
optimization  driver  (QNMDIF)  was  carried  out.  As  could  be  anticipated,  the 
resulting  combined  code  is  quite  large,  i.e.  approximately  17,000  lines  of 
Fortran  source.  Consequently,  we  have  decided  that  in  order  to  expedite  the 
initial  testing  of  the  approximation  method  in  this  combined  procedure  that  the 
segment  of  the  code  containing  the  approximation  method  together  with  all  the 
geometry  and  other  routines  that  are  required  to  enable  prediction  of  the 
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approximation-predicted  objective  function  be  separately  modularized.  This 
was  done  in  such  a  fashion  as  to  allow  separate  testing  of  the  ^proximation 
methods'  performance  in  the  solution  space  typical  of  the  optimization  design 
environment,  and  then  to  later  allow  direct  coupling  with  the  QNMDIF/TWING 
optimization  procedure  when  the  approximation  method  testing  is  complete. 

3.3.2  Definition  of  3-D  Turbomachinery  Design  Optimization  Problem  for 
Testbed  Case  Studies  With  Approximation  Method 

After  reviewing  the  various  classes  of  optimization  problems  that  can  be  directly 
implemented  with  the  QNMDIF/TWING  optimization  procedure  as  it  is  presently 
constituted,  a  particular  class  of  aerodynamic  optimization  design  problems  was 
selected  that  is  sufficiently  general  in  its  overall  design  aspects  to  be 
representative  of  practical  design  problems  relevant  to  turbomachinery 
applications.  The  selected  problem  class  relates  to  the  optimized  design  at 
highly  supercritical  flow  conditions  of  a  3-D  wing  or  isolated  blade  geometry  so 
as  to  minimize  its  total  drag-to-lift  ratio.  This  choice  provides  a  challenging 
design  problem  to  the  optimization  method  as  well  as  the  approximation 
method  because  of  both  the  overall  sensitivity  of  the  3-D  flow  to  small  changes 
in  geometry,  as  well  as  the  particular  choice  of  objective  function.  Several 
investigators  have  found  (Refs.  23-25)  that  although  the  selection  of  drag  or 
drag/lift  ratio  as  the  objective  function  is  the  most  relevant  choice  to  the  design 
process,  drag  was  too  sensitive  a  quantity  for  the  optimization  and  inviscid  flow 
solver  codes  employed  at  that  time  to  produce  acceptable  results.  In  order  to 
achieve  some  measure  of  optimization,  the  less  satisfactory  choice  of 
employing  a  target  pressure  distribution  as  the  objective  had  to  be  made.  From 
a  physical  standpoint,  this  is  not  nearly  as  desirable  as  using  overall  drag  as  an 
objective  since  apriori  knowledge  of  an  appropriate  target  pressure  distribution 
is  not  conveniently  available  for  most  problems.  However,  it  has  been  found 
(Ref.  8)  that  by  using  the  new  quasi-Newton  optimization  driver  together  with  an 
accurate  full  potential  flow  solver  and  a  higher  3-D  flow  solution  convergence 
tolerance  the  employment  of  inviscid  drag  as  the  objective  function  is  practical, 
and  that  good  optimization  results  can  be  achieved  on  that  basis. 

Consequently,  based  upon  these  demonstrated  results  the  initial  testbed 
optimization  case  study  for  the  3-D  approximation  method  was  selected  to  be 
the  optimized  design  of  a  3-D  wing  or  isolated  blade  operating  at  highly 
supercritical  flow  conditions  using  the  drag/lift  ratio  as  the  objective  function  with 
approximately  10  design  variables  related  to  wing  surface  geometry.  This  is  a 
challenging  design  environment  for  the  approximation  method  since  the  3-D 
flow  solutions  required  here  are  extremely  nonlinear,  and  the  simultaneous 
multiple  design  variable  parameter  changes  required  by  the  optimizer  as  it 
proceeds  through  the  design  variable  solution  space  usually  exhibit  a  strong 
nonlinear  coupling  between  the  solution  states.  However,  the  overall  design 
problem  represented  by  this  selection  is  so  fundamental  and  directly  relatable 
to  practical  aerodynamic  design  that  its  accurate  solution  by  the  approximation 
method  would  be  a  major  contribution  to  achieving  the  general  transition  of 
turbomachinery  analysis  codes  to  design  mode  utilization. 
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The  specific  design  strategy  selected  is  based  on  that  employed  in  Ref.  8  since 
it  appears  to  be  one  of  the  most  efficient  optimization  strat^ies  reported  to  date. 
The  design  problem  involves  the  geometric  optimization  of  the  surface  of  a  3-D 
wing  or  isolated  turbomachinery  blade  using  a  relatively  small  number  of 
design  variables  to  modify  the  surface  geometry.  The  overall  problem  is 
outlined  in  Figure  1 5.  In  contrast  to  previous  optimization  studies  which  have 
employed  various  analytic  shape  functions  distributed  along  the  chord  of  the 
profile  (Refs.  20-22)  to  alter  the  baseline  profile,  the  technique  suggested  in  Ref. 
8  is  to  simply  assign  as  design  variables  the  vertical  displacement  of  several 
pre-assigned  points  on  the  wing  surface.  Then,  as  shown  in  Figure  15,  by 
holding  invariant  the  surface  geometry  along  that  portion  of  the  profile  surface 
we  do  not  wish  to  change,  the  surface  geometry  along  the  remainder  of  the 
profile  is  altered  by  employing  a  spline-fitting  routine  to  redefine  that  portion  of 
the  surface  as  the  pre-assigned  movable  points  displace  vertically  from  their 
baseline  locations  in  response  to  the  optimization  process.  The  advantage  of 
this  strategy  is  that  a  small  number  of  movable  points  is  sufficient  to  redesign  a 
relatively  large  segment  of  a  3-D  geometric  surface.  For  instance,  for  the 
example  shown  in  Figure  1 5,  at  each  of  3  spanwise  locations  on  the  wing  three 
chordwise  points  located  at  approximately  x/c  =  {0.2,0.4,0.6}  were  found  to  be 
sufficient  to  recontour  the  upper  surface  segment  of  the  wing  along  the 
chordwise  band  from  x/c  =  {0.08  to  0.75}  and  to  enable  quite  satisfactory 
optimization  results  to  be  obtained.  We  will  discuss  these  results  shortly.  In 
contrast,  employing  analytic  shape  functions  distributed  at  various  locations 
along  the  profile  chord  would  typically  require  at  least  double  the  number  of 
design  variables. 

The  initial  three-dimensional  optimization  case  study  chosen  for  the 
approximation  method  is  then  summarized  as  follows:  using  drag/lift  as  the 
objective,  the  planform  geometry  as  illustrated  in  Figure  15,  and  a  constant 
profile  geometry  across  the  span  given  by  the  NACA  64  A212  airfoil,  the  upper 
surface  of  the  wing  was  modified  along  the  chordwise  band  from  approximately 
x/c  =  {0.08  to  0.75}  using  a  total  of  9  design  variables  associated  with  the 
movable  spline-supported  points  as  indicated  Figure  15  and  located  at  the  3 
spanwise  locations  y/s  =  {0.25,0.50,.075}  and  chordwise  locations  x/c  » 
{0.2,0.4,0.6}. 


3.3.3  Results  of  Combined  Approximation  Method  and  Optimization  Procedure 
for  3-D  Design  Optimization  Case  Studies 

In  order  to  establish  the  benchmark  for  evaluating  the  approximation  methods 
performance,  we  have  run  the  ONMDIF/TWING  optimization  procedure  without 
employing  the  approximation  method.  This  particular  optimization  problem  was 
previously  studied  in  Ref.  8  so  that  a  direct  comparative  result  is  available. 
Figure  16  provides  the  original  spanwise  pressure  distributions  of  the  baseline 
configuration,  and  provides  in  the  lower  left  of  the  figure  the  nominal  values  of 
the  nine  design  variables  for  the  baseline  configuration  multiplied  by  100,  i.e. 
values  of  the  normalized  Z  coordinate  of  the  movable  spline-fitted  points  100-{ 
(Z/c)^,  (Z/c)2^  ...}.  Figures  17,  18,  19,  and  20  display  the  optimized  wing 
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spanwise  pressure  distributions  after  the  7th,  8th.  9th,  and  10th  optimization 
cycles,  respectively.  The  10  cycles  required  approximately  12,000  secs,  of 
CRAY  XMP  CPU  time.  In  our  calculations  we  have  employed  a  finer  grid 
(149,31,26)  than  that  (89,25,18)  previously  used  since  we  wished  to  resolve  the 
upper  surface  shock  more  accurately.  Figure  21  provides  an  interesting 
synopsis  of  the  behavior  of  the  quasi-Newton  optimization  procedure  for  this 
design  problem.  That  figure  provides  the  behavior  history  of  the  objective 
function  each  time  it  is  evaluated  from  the  3-D  flow  solution  as  called  by  the 
optimization  driver.  As  can  be  seen,  essentially  the  first  30  evaluations  and 
corresponding  3-D  flow  solutions  are  needed  to  set  up  the  Hessian  matrix  of 
^  mixed  second  derivatives  of  the  objective  function  with  respect  to  the  9  design 

variables.  Once  that  is  completed,  there  is  a  large  drop  in  objective  function  as 
the  method  proceeds  to  its  first  design  point.  After  that,  only  gradient 
evaluations  are  needed,  i.e.  essentially  one  3-D  solution  for  each  design 
variable,  to  take  the  method  through  the  succesive  optimization  cycles.  The 
rapid  and  continuing  convergence  of  the  procedure  is  evident,  in  contrast  to  the 
slow  and  often  halting  convergence  behavior  of  the  conjugate  gradient  method. 

As  the  first  test  of  the  3-D  approximation  method's  ability  to  predict  3-D  surface 
pressure  distributions  in  the  multiple  design  variable  solution  space  associated 
with  the  above  optimization  problem  we  have  employed  the  approximation 
method  to  predict  the  surface  pressure  distributions  for  the  design  variable 
values  associated  with  the  design  point  at  the  end  of  each  of  the  last  4 
optimization  cycles  run  on  the  above  design  problem.  These  are  the  pressure 
distributions  previously  shown  in  Figures  17-20.  Figures  22,  23,  24,  and  25, 
respectively,  provide  comparisons  of  the  approximation  predicted  surface 
pressures  at  the  wing  root  chord  with  the  exact  3-D  TWING  flow  solutions  at 
those  same  4  design  points  illustrated  in  Figures  17-20.  With  the  exception  of 
several  points  in  the  near  vicinity  of  the  upper  surface  shock  wave,  the 
comparisons  are  quite  reasonable.  In  the  numerical  determination  of  the  3-D 
flow  solution  by  TWING,  the  computational  grid  used  employs  25  spanwise 
locations  on  the  wing.  At  all  of  these  25  spanwise  stations,  an  approximation 
prediction  of  surface  pressure  is  required  in  order  to  accurately  calculate  the 
approximation-predicted  drag  for  comparison  with  the  exact  result.  For  this  test 
case,  however,  at  approximately  the  sixth  spanwise  location  we  found  the  initial 
approximation  method  predictions  to  break  down  due  to  shock  wave  invariant 
point  crossings.  This  is  not  entirely  surprising  considering  the  topologically 
complex  pressure  distribution  presented  by  the  current  optimization  problem. 
Note  that  the  approximation  method  is  trying  to  track  both  the  chordwise  and 
spanwise  movement  of  the  shock  wave  due  to  the  simultaneous  change  of  9 
•  design  variables,  all  occurring  in  a  highly  nonlinear,  coupled  flow  environment. 

What  is  needed  to  remedy  this  occurrence  in  the  approximation  method  is  a 
more  generalized  and  robust  determination  of  the  3-D  invariant  point  locations 
that  are  required  by  the  approximation  method.  We  discuss  its  development 
below. 

Another  issue  of  note  is  the  strategy  for  selecting  the  calibration  solution  matrix 
required  as  input  to  the  approximation  method.  For  our  first  choice,  we  simply 
selected  a  subset  of  the  3-D  solutions  determined  by  the  QNMDIF  optimization 
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method  for  the  Hessian  matrix  calculation.  This  strategy  has  the  appeal  that  no 
additional  3-D  flow  solutions  are  needed  to  be  calculated.  What  we  found, 
however,  was  that  the  design  variable  stepsizes  initially  selected  by  the 
optimizer  are  too  small  and  do  not  cover  a  sufficient  volume  of  the  design 
variable  solution  space  that  will  ultimately  be  searched  by  the  optimizer  as  it 
proceeds  to  the  optimum  design  point.  Consequently,  for  the  approximation 
predicted  results  presented  in  Figures  22  to  25  the  calibration  solution  matrix 
was  hand  picked  so  as  to  provide  a  reasonable  coverage  of  the  flow  solution 
space  from  the  baseline  solution  shown  in  Figure  16  to  the  final  design  point 
solution  shown  in  Figure  20.  Of  course,  apriori  knowledge  of  the  appropriate 
range  for  design  variable  movement  is  generally  not  available.  Consequently, 
the  most  effective  strategy  for  determining  the  calibration  solution  matrix,  and 
perhaps  even  more  generally,  the  best  way  to  employ  the  approximation 
method,  i.e.  from  the  very  start  of  the  optimization  problem  so  that  even  the 
Hessian  matrix  is  determined  using  approximation-predicted  3-D  solutions,  or 
immediately  after  the  Hessian  matrix  is  determined,  or  at  the  end  of  the  first 
optimization  cycle,  or  some  combination  or  variation  of  these  different 
strategies,  requires  investigation. 

In  working  with  the  combined  QNMDIF/TWING/APPROX  pilot  code.  It  was 
immediately  found  to  be  necessary  to  separate  and  mcKlularize  the  3-D 
geometry  and  grid  generator  routines  from  the  combined  pilot  code  so  as  to 
allow  direct  interaction  with  either  full  3-D  TWING  solutions  or  with 
approximation-predicted  solutions.  The  reason  for  separating  and  modularizing 
the  3'D  wing  geometry  and  finite  difference  grid  generator  from  the  combined 
QNMDIF/TWING/APPROX  pilot  code  was  to  permit  the  extensive  testing  of  the 
approximation  method  to  be  carried  out  without  continually  being  burdened  by 
the  overhead  of  the  entire  QNMDIF/TWING/APPROX  program  combination. 
Carrying  this  overhead  is  not  necessary  for  testing  the  approximation  method, 
and  modularization  of  these  routines  adds  considerably  to  the  computational 
efficiency.  Note  that  just  as  with  the  full  3-D  TWING  solutions,  in  order  to 
compute  the  approximation-predicted  lift  and  drag,  it  is  necessary  to 
redetermine  the  wing  geometry  for  each  variation  in  design  variables  so  that 
pressure  integrations  can  be  carried  out  over  the  newly-designed  wing  surface. 
In  determining  full  3-D  TWING  solutions,  the  TWING  grid  generator 
redetermines  the  chordwise  distributions  of  wing  surface  points  that  are  initially 
determined  by  the  wing  geometry  generator  and  then  input  to  it.  We  also 
elected  to  employ  the  TWING  grid  generator  to  determine  the  points  at  which 
the  approximation-predicted  pressures  were  to  be  determined.  This  insures 
that  for  the  approximation-predicted  results  a  smooth  distribution  of  chordwise 
points  will  be  achieved  at  each  spanwise  location.  Additionally,  it  further 
guarantees  that  unnecessary  inaccuracies  associated  with  solution 
interpolations  are  avoided  when  comparing  approximation-predicted  and  exact 
TWING  results. 

After  the  separation  and  modularization  of  those  geometry  and  grid  generator 
routines,  verification  calculations  of  the  lift  and  drag  using  the  approximation 
method  were  performed  in  order  to  verify  the  modularization.  This  testing  was 
accomplished  by  employing  the  modularized  geometry  and  grid  generator 
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routines  together  with  the  approximation  method  to  reproduce  the  results 
associated  with  the  3-D  TWING  base  and  calibration  solutions.  That  was 
carried  out  by  setting  the  design  variable  values  to  be  those  of  the  base  and 
calibration  solutions.  Hence,  the  approximation  method  then  simply  repredicts 
the  base  and  calibration  pressure  distributions,  if  the  geometry  and  grid 
generator  routines  then  correctly  determine  the  wing  geometry  information 
required  to  calculate  the  lift  and  drag  at  each  spanwise  location,  the  results  so 
calculated  should  compare  identically  with  those  detennined  from  the  exact  3-D 
TWING  calculations.  It  was  verified  that  the  predicted  sectional  lift  and  drag 
predicted  by  the  modularized  code  exactly  reproduces  that  determined  by  the 
full  TWING  code. 

Examples  of  the  output  of  the  predicted  sectional  lift,  drag,  and  drag/lift  ratio  as 
determined  from  the  modularized  code  are  provided  in  Figures  26-29.  For 
example,  illustrated  in  Figure  26a.  are  the  spanwise  distributions  at  each  of  the 
25  spanwise  locations  along  the  wing  of  the  sectional  lift,  drag,  and  drag/lift  ratio 
for  the  optimized  wing  surface  pressure  distribution  shown  in  Figure  26b.  Note 
that  this  pressure  distribution  is  the  result  that  was  obtained  after  the  7th 
optimization  cycle  for  the  benchmark  case  study. 

Additional  successfully-predicted  results  from  the  modularized  code  for  the 
sectional  lift  and  drag  are  provided  in  Figures  27a.-b.  to  29a.-b.  Those  results 
are  for  the  optimized  pressure  distributions  after  the  8th.  9th.  and  10th 
optimization  cycles,  resp^ively,  for  the  benchmark  case  study.  Consequently, 
the  results  illustrated  in  Figures  26-29  represent  both  the  target  pressure 
distributions  and  the  sectional  lift  and  drag  results  that  the  approximation 
method  will  attempt  to  predict  for  the  testbed  optimization  study. 


3.4  Extended  Invariant  Point  Relocation  Procedure  and  Postprocessing  of 
Approximation  Results 

An  improved  methodology  for  redefining  the  3-D  chordwise  and  spanwise 
invariant  point  locations  associated  with  shock  waves  that  occur  in  the  base  and 
calibration  solutions  was  developed.  This  was  found  to  be  necessary  as  a 
result  of  testing  carried  out  on  the  benchmark  optimization  case  study.  This 
involved  determination  of  a  procedure  for  relocating  the  shock  wave  invariant 
points  in  the  spanwise  direction  that  provides  a  continuous  and  smooth 
variation  of  both  the  invariant  point  locations  as  well  as  the  associated  values  of 
the  pressure  at  those  points.  An  illustration  of  the  typical  improvement  provided 
by  the  relocation  procedure  is  provided  in  Figures  ^a.-b.  to  32a.-b.  In  Figure 
30a.,  results  are  given  for  the  unmodified  base  TWING  solution  employed  in  the 
benchmark  case  study.  The  plot  on  the  left  of  that  figure  displays  the  spanwise 
locations  of  the  invariant  points  associated  with  the  upper  surface  shock,  while 
the  plot  on  the  right  shows  the  associated  Cp  values  of  those  same  points.  In 
both  of  the  plots,  the  triangular  symbols  represent  the  Invariant  point  locations 
associated  with  the  beginning  of  the  shock  wave,  while  the  square  and  circular 
symbols  represent,  respectively,  the  correponding  sonic  point  and  shock 
termination  locations.  Note  the  nonsmoothness  evidenced  by  several  of  the 
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points,  particularly  as  shown  in  the  Cp  values.  Figure  30b.  illustrates  the  same 
plots  after  the  relocation  procedure  has  been  applied.  The  smooth  variation  in 
both  location  and  Cp  values  is  evident.  More  remarkable  results  of  the 
procedure  are  provided  in  Figure  31a.'b.  That  figure  illustrates  analogous 
results  for  one  of  the  calibration  solutions  in  which  the  original  invariant  point 
location  procedure  could  not  locate  a  distinguishable  invariant  point  marking 
the  beginning  of  the  shock  wave  at  several  spanwise  locations  near  the  root 
chord,  as  can  seen  in  Figure  31a.  However,  as  shown  in  Figure  31b.,  the 
relocation  procedure  successfully  provides  smoothly-varying  invariant  points 
marking  the  beginning  of  the  sho^  wave  at  those  lotions  for  this  calibration 
solution.  The  final  result  provided  in  Figures  32a.’b.  illustrates  how  irregular  the 
locations  of  the  invariant  locations  associated  with  certain  shock  topologies  can 
become.  The  plot  shown  in  Figure  32a.  exhibits  a  gap  occurring  at 
approximately  mid-span  in  the  identification  of  the  invariant  points  associated 
with  the  pre-shock  and  sonic  locations  in  the  pressure  distribution  for  this 
calibration  solution.  This  phenomena  has  apparently  occurred  since  this 
particular  pressure  distribution  does  not  recover  through  sonic  velocity  over  that 
spanwise  location.  Although  while  this  phenomena  may  occur  physically,  such 
a  topology  change  will  prevent  a  prediction  by  the  approximation  method  if  not 
corrected.  The  relocation  procedure  corrects  this  problem  by  beating  a  smooth 
distribution  of  invariant  points  across  that  spanwise  gap  as  shown  in  Figure 
32b. 

Finally,  as  a  result  of  all  the  above  development  of  the  3-D  approximation 
method.  Figures  33a.-y.  display  a  complete  prediction  of  a  3-D  pressure 
distribution  by  the  approximation  method  as  it  is  presently  constituted. 
Illustrated  in  these  figures  are  comparisons  of  the  approximation  method  and 
the  exact  TWING  results  for  the  chordwise  pressure  distributions  at  all  of  the  25 
spanwise  stations  across  the  wing  at  whbh  the  grid  generator  has  located  the 
computational  grid  for  the  3-D  pressure  distribution  associated  with  the  7th 
optimization  cycle  of  the  benchmark  ce^e  study  (see  Figure  26b.).  Also  shown 
on  each  of  those  figures  for  comparative  purposes  are  the  base  and  9 
calibration  pressure  distributions  which  are  employed  by  the  approximation 
method  to  predict  the  pressure  distribution  results  shown.  The  corresponding 
prediction  of  the  spanwise  sectional  lift,  drag,  and  drag/lift  ratio  are  given  in 
Figure  34a.  and  compared  with  the  exact  result  in  Figure  34b.  taken  from  Figure 
26.  As  can  be  seen  in  Figure  34.  with  the  exception  of  spanwise  bcations  15. 
19,  and  25.  quite  reasonable  predictions  of  these  integrated  spanwise 
quantities  are  obtained  from  the  approximation-predicted  pressures.  However, 
in  order  to  insure  a  consistent  and  smooth  spanwise  variation  of  all  the 
integrated  quantities,  development  of  a  post-processing  procedure  to  correct 
and  rebcate  the  final  ^proximation-predicted  invariant  point  bcations  in  the 
spanwise  direction  is  still  necessary  and  is  discussed  below. 

Considering  Figure  34  and  setting  aside  for  the  moment  the  fact  that  the  total 
integrated  lift,  drag,  and  drag/lift  ratio  of  the  approximation-predicted  results 
indicated  in  Figure  34a.  are  somewhat  in  error  compared  with  the  exact  result 
provided  in  Figure  34b.,  the  overall  approximation  prediction  for  the  spanwise 
variation  of  sectional  quantities  appears  quite  reasonable.  In  our  initial 
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investigation  of  this  problem,  it  appeared  that  the  major  source  of  error 
preventing  accurate  prediction  by  the  approximation  method  of  the  total 
integrated  forces  lie  with  the  several  spanwise  locations,  such  as  locations  15, 
19,  and  25  shown  in  Figure  34a..  in  which  clearly  inconsistent  behavior  is 
observed  in  the  approximation  prediction  of  the  sectional  quantities. 
Implementation  of  an  additional  post-processing  invariant-point  smoothing 
procedure  for  correcting  the  final  approximation-predicted  invariant  point 
locations  in  the  spanwise  direction  appeared  straightforward  and  was  carried 
out.  Comparison  of  the  post-processed  invariant  point  spanwise  pressure 
distributions,  which  are  not  shown  here,  with  those  in  Figure  34a.  exhibited 
some  variation  from  the  nonpost-processed  distributions  but  resulted  in  no 
dramatic  changes.  Figure  35  provides  the  new  approximation  results  based  on 
post-processed  spanwise  surface  pressure  distributions  to  be  compared  with 
the  previous  approximation  and  exact  results  shown  in  Figures  34a.  and  34b., 
respectively.  This  new  approximation  result  exhibits  somewhat  stronger 
discontinuous  behavior  at  several  spanwise  locations  than  other  results  that  we 
have  obtained  when  applying  slightly  different  spanwise  smoothing  criteria  in 
the  invariant  point  location  procedure,  and  it  illustrates  the  basic  sensitivity 
problem  that  we've  encountered  with  these  multi-parameter  3-D  solution 
predictions.  That  is.  that  while  the  invariant  point  relocation  procedures 
developed  and  employed  here  with  the  3-D  approximation  method  can  and  do 
insure  continuous  variations  of  invariant  point  locations  in  the  spanwise 
direction,  a  continuous  variation  in  the  sectional  integrated  quantities  obtained 
from  those  pressure  distributions,  such  as  lift  and  drag,  is  not  similarly  insured. 
The  reason  for  this  is  due  to  the  strong  tendency  for  extremely  rapid  local 
topology  change  in  the  basic  surface  pressure  distribution  for  these  strongiy- 
nonlinear  flows.  These  rapid  topology  changes  in  surface  pressure  can  have 
the  result  that  a  small  error  in  the  spanwise  invariant  point  location  distribution 
can  cause  the  strength  and  location  of  the  shock  at  one  or  more  spanwise 
locations  to  be  sufficiently  in  error  to  cause  a  significant  discontinuous  variation 
in  the  spanwise  variation  of  the  integrated  quantities. 

The  most  direct  way  of  correcting  for  this  problem  is  to  further  enhance  the 
current  invariant  point  relocation  procedure  both  for  redefining  the  invariant 
point  locations  in  the  base  and  calibration  solutions  as  well  as  post-processing 
the  final  approximation-predicted  invariant  point  distributions.  An  alternative 
would  be  the  possibility  of  working  with  another  basic  flow  property  rather  than 
such  a  strongly  sensitive  quantity  as  surface  pressure.  An  interesting  idea 
which  comes  to  mind  is  that  instead  of  employing  the  pressure  as  our 
fundamental  dependent  quantity  which  is  discontinuous  through  a  shock,  we 
employ  instead  the  mass  flux  or  some  similar  quantity  that  remains  continuous 
through  a  shock,  and  then  subsequently  determine  the  flow  properties  that  are 
of  interest  from  these  quantities. 

Recall  that  the  primary  function  of  the  invariant  point  location  procedure  is  to 
provide  a  continuous  spanwise  variation  of  both  the  invariant  point  locations  as 
well  as  the  associated  values  of  the  pressure  at  those  points.  The  relocation 
procedure  corrects  descrepancies  due  to  numerical  'jitter*  in  the  vicinity  of  high- 
gradient  regions  such  as  shock  waves  in  the  base  and  calibration  solutions. 
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This  jitter  if  uncorrected  results  in  significant  degradation  of  the  unit  perturbation 
distributions  which  are  determined  by  differencing  in  strained  coordinates  the 
various  calibration  solution  distributions  and  the  base  solution  distribution,  in 
addition  to  employing  the  invariant  point  relocation  procedure  on  the  base  and 
calibration  solutions,  which  can  be  thought  of  as  preprocessing  those  solutions 
in  preparation  for  use  by  the  approximation  solution  procedure,  as  a  result  of 
additional  numerical  exf^rimentation  with  the  approximation  method  we  have 
also  found  it  necessary  to  again  apply  an  invariant  point  relocation  procedure  to 
the  approximation-predicted  invariant  points.  This  in  effect  postprocesses  the 
approximation  solution.  We  have  found  this  to  be  necessary  due  to  the  fact  that 
the  approximation-predicted  location  of  the  invariant  points  often  contains 
numerical  jitter  as  well,  particularly  for  multi-parameter  perturbations  involving 
four  or  more  simultaneous  perturbations.  This  jitter  often  results  in 
nonuniformities  in  the  spanwise  distribution  of  the  invariant  points.  Although  the 
spanwise  continuity  of  the  predicted  pressure  distributions  can  appear 
reasonable,  integration  of  these  distributions  across  the  span  direction  often 
displays  nonuniform  variation  in  the  integrated  sectional  quantities  such  as  lift 
and  drag. 

In  order  to  investigate  the  underlying  cause  and  determine  the  appropriate 
method  to  correct  the  observation  of  the  descrepancies  in  the  spanwise 
variation  of  the  approximation-predicted  sectional  lift  and  drag,  we  have 
performed  a  detailed  evaluation  of  the  performance  of  the  approximation 
method  involving  one,  two,  three,  and  nine  design  parameters  when  employing 
combined  preprocessing  of  the  invariant  points  in  the  base  and  calibration 
solutions  plus  postprocessing  of  the  invariant  points  in  the  3-D  approximation- 
predicted  solution.  As  a  result  of  this  numerical  testing,  we  have  made  the 
following  observations.  Application  of  the  invariant  point  relocation  procedure 
to  smooth  out  the  numerical  jitter  in  the  invariant  point  locations  in  the  original 
base  and  calibration  solutions  is  very  effective  and  results  in  accurate 
approximation  predictions  of  the  new  invariant  point  locations  for  situations 
involving  a  small  number  of  simultaneous  multiple  parameter  perturbations.  We 
have  found  that  for  the  cases  we  have  consider^,  a  two  or  three-parameter 
perturbation  problem  can  usually  be  accurately  modeled  by  the  approximation 
procedure.  The  tested  problems  involved  strongly  supercritical  flows  with  a 
strong,  single  shock  wave  across  the  upper  wing  surface.  However,  when  the 
number  of  multiple  parameters  to  be  simultaneously  altered  increases  beyond 
three  then  what  seems  inevitably  to  occur  in  the  approximation  prediction  is  that 
some  small  numerical  jitter  in  the  invariant  point  locations,  steming  from  either 
imperfections  in  the  originai  base  and  calibration  solution  locations  or  from  the 
preprocessing  procedure  that  may  slightly  misalign  one  or  more  invariant 
points,  results  in  a  large  discrepan<7  in  the  predicted  resuits  at  one  or  more  of 
the  spanwise  locations.  In  order  to  treat  this  problem,  we  have  proceeded  to 
develop  and  test  a  procedure  to  relocate  the  invariant  points  in  the  predicted 
solution  prior  to  using  the  approximation-predicted  surface  pressure 
distributions  to  determine  integrated  sectional  quantities.  As  a  result  of  this,  we 
have  been  able  to  demonstrate  that  the  employment  of  a  post-processing 
procedure  to  relocate  and  realign  the  various  invariant  point  distributions  in  the 
approximation-predicted  results  can  resolve  essentially  all  of  the  discrepancies 
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in  the  predicted  results  again  when  the  number  of  varied  parameters  is  small. 
However,  when  the  number  of  multiple  parameters  increases  to  a  larger 
number  typically  required  in  a  realistic  design  problem,  a  number  which  would 
usually  be  10  or  more,  the  current  smoothing  procedures  cannot  guarantee 
consistent  results.  The  basic  requirement  on  the  methods  being  developed  in 
this  study  is  to  track  and  accurately  replicate  the  motion  of  invariant  points 
characterizing  high-gradient,  discontinuous  regions  of  the  flow  field.  The 
prediction  capability  of  these  methods  is  critically  dependent  upon  the 
goodness  of  the  basic  numerical  solutions  for  the  base  and  calibration  flow 
fields  that  are  input  to  them.  For  classes  of  highly  sensitive  flows  as  we  have 
been  studying,  which  are  inevitably  the  ones  of  most  interest  for  practical  design 
applications  since  they  usually  represent  situations  that  can  benefit  the  greatest 
from  design  refinements,  small  discrepancies  in  the  basic  numerical  flow 
solutions  that  are  input  to  the  approximation  procedure  can  lead  to  large 
discrepancies  in  the  predicted  pressure  distributions  in  these  high-gradient 
regions.  Since  some  numerical  jitter  in  these  regions  in  the  base  and 
calibration  numerical  solutions  is  unavoidable,  the  resultant  accuracy  and 
range  of  validity  of  the  approximation  method  predictions  will  be  significantly 
reduced  in  these  applications  when  the  basic  input  to  the  approximation 
method  consists  of  highly-sensitive  surface  pressure  distributions. 

An  alternative  method  to  alleviate  this  basic  restriction  on  the  accuracy  of  the 
approximation  method,  as  mentioned  briefly  above,  is  to  employ  as  input  to  the 
approximation  method  flow  properties  other  than  those  on  the  wing  or  blade 
surface  where  the  gradients  are  naturally  the  highest.  Since  most  practical 
design/optimization  problems  in  both  aerodynamics  and  turbomachinery  are 
focused  on  seeking  an  improvement  of  an  integrated  flow  quantity,  such  as  drag 
or  loss  coefficient,  working  with  surface  properties  is  not  always  necessary.  For 
example,  in  order  to  determine  the  forces  acting  upon  an  aerodynamic 
configuration  one  can  employ  control  surfaces  removed  from  the  actual 
aerodynamic  surface  and  located  at  any  convenient  location  in  the  flow  field 
where  information  is  available.  Because  numerical  flow  field  solution  methods 
provide  information  everywhere  in  the  flow  field,  that  information  is  already 
available,  in  the  present  study,  as  well  as  in  all  of  the  previous  work  involving 
the  approximation  method,  we  have  employed  surface  pressure  as  the  basic 
dependent  variable  since  surface  pressure  is  both  a  fundamental  aerodynamic 
quantity  and  surface  pressure  distributions  were  found  to  be  very  convenient  to 
use  in  the  comparative  testing  involved  In  the  preliminary  development  of  the 
approximation  method.  However,  there  is  no  fundamental  need  or  restriction  to 
employing  surface  quantities  with  the  approximation  method.  Any  convenient 
distributions  will  serve  equally  as  well.  For  our  purposes,  it  appears  just  as 
suitable  to  employ  distributions  of  flow  properties  at  the  outer  edge  of  the 
computational  grid  where  the  flow  contains  significantly  less  high-gradient 
features  and  is  more  of  a  small-disturbance  flow.  A  strict  small  disturbance  flow 
Is  not  necessary  for  accurate  use  with  the  approximation  method.  The  far-fiekJ 
flow  can  contain  one  or  several  shocks  passing  through  the  outer  boundary  and 
the  approximation  method  will  work  accurately  without  restriction.  The  essential 
point  is  that  the  features  of  the  flow  field  at  points  removed  from  the  actual  blade 
surface  will  be  less  severe  and  subsequently  less  susceptable  to  numerical 
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jitter  than  for  example  at  the  foot  of  a  strong  shock  wave  on  the  surface.  In  fact, 
employment  of  outer  flow  field  information  to  predict  integrated  quantities  has 
been  recognized  in  the  past  and  previously  employed.  Henne  at^  Hicks  (Ref. 
29)  have  developed  and  employed  a  contour  integration  procedure  applied  to 
the  outer  computational  grid  boundary  of  the  3-D  transonic  flow  code  FL022NM 
to  evaluate  the  lift  and  drag  of  wings  and  wing/bodies  at  transonic  speeds.  This 
was  done  in  order  to  avoid  working  with  surface  pressures  to  determine  the 
drag  which  was  found  to  be  in  error  due  to  small  numerical  inaccuracies  in  the 
predicted  surface  pressures.  The  result  was  a  consistent  and  more  accurate 
prediction  of  drag  throughout  the  transonic  regime.  We  suggest  that  the 
alternative  of  employing  flow  field  rather  than  surface  properties  as  the  basic 
input  to  the  approximation  method  be  examined. 

Before  making  the  final  decision  on  whether  to  alter  the  mode  of  operation  of 
the  approximation  method  from  employing  surface  pressure  distributions  as  the 
basic  input  information  to  using  flow  field  properties,  we  have  carefully  reviewed 
all  of  our  previous  approximation  method  results  involving  single  and  multiple- 
parameter  simultaneous  changes  to  detemiine  the  precise  sources  that  cause 
the  discrepancies  currently  observed  in  the  predicted  results  and  to  insure  that 
no  simple  mechanism  for  easily  correcting  these  discrepancies  has  been 
overlooked.  We  have  reviewed  the  following  four  major  steps  in  the 
approximation  procedure:  (1)  the  preprocessing  procedure  for  identifying  and 
relocating  the  invariant  point  locations  in  the  base  and  calibration  solutions  in 
order  to  provide  a  smooth  variation  of  these  invariant  points  in  the  spanwise 
direction,  (2)  the  procedure  for  determining  the  approximation-predicted 
solution  based  upon  the  preprocessed  base/calibration  solutions,  in  particular 
the  prediction  of  the  new  invariant  point  locations,  (3)  the  post-processing 
procedure  for  relocating  the  predicted  invariant  point  locations  in  the  spanwise 
direction,  (4)  the  integration  of  the  post-processed  approximation-predicted 
solutions  to  determine  integrated  sectional  quantities  in  the  spanwise  direction. 

For  the  class  of  flow  solutions  characterizing  the  benchmark  case  study,  it 
appears  that,  even  after  careful  post-processing,  integrated  sectional  results  of 
the  approximation-predicted  solutions  tend  to  exhibit  noncontinuous  spanwise 
behavior  due  to  invariant  point  misalignment.  Although  the  benchmark  test 
case  study  involves  topologically  complex  flows  in  the  sense  that  strongly 
supercritical  flows  are  involved,  we  believe  it  to  be  a  good  test  of  the 
approximation  method  if  the  method  is  to  achieve  a  reasonably  general 
applicability  to  turbomachinery  flows  where  flow  topologies  are  at  least  as 
complex  as  those  considered  here. 

As  a  result  of  our  review,  we  believe  we  have  now  identified  the  multiple 
sources  causing  the  limitation  in  the  application  of  the  approximation  method  in 
three-dimensional  studies.  The  primary  factors  are  a  combination  of  numerical 
inaccuracies  in  the  base/calibration  solutions  near  the  invariant  point  locations 
associated  with  the  foot  of  the  shock  together  with  a  geometrical  compounding 
of  these  errors  when  the  approximation-predicted  invariant  point  locations  are 
dependent  upon  the  simultaneous  variation  of  a  large  number  of  multiple 
parameters.  The  fundamental  action  of  the  approximation  method  is  to  move 
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the  invariant  points  lineariy  in  multiply-strained  coordinate  space  with  respect  to 
each  perturbation  parameter.  The  predicted  results  from  such  a  formulation 
have  been  found  to  be  extremely  accurate  in  2-D  situations  where  the  shock 
structure  was  relatively  simple  and  well-resolved  by  the  numencal  solution 
procedure  so  that  the  associated  invariant  points  were  also  well  defined.  For 
three-dimensional  situations,  however,  the  shock  structure  on  or  near  the 
surface  of  the  3-D  configuration  is  not  so  well  defined  by  the  flow  solvers. 
Consequently,  the  definition  of  the  invariant  point  locations  in  the 
base/caiibration  solutions  is  correspondingly  lacking.  Attempting  to  correct 
these  locations  by  employing  general  smoothing  criteria  can  improve  the 
situation  at  some  locations  but  inevitably  increases  the  error  at  other  locations. 
Finally,  when  the  approximation-predicted  invariant  point  locations  are 
determined  as  a  sum  of  all  the  individual  contributions  from  each  varied 
parameter,  the  possibility  of  a  significant  spurious  result  occurring  at  even  one 
location  is  high.  For  example,  in  a  problem  involving  10  design  variables,  if  the 
approximation  prediction  accurately  accounts  for  the  invariant  point  motions 
and  resultant  flow  topology  changes  of  9  of  the  design  variables,  but  misses  the 
prediction  for  the  10th  at  one  spanwise  location,  the  resultant  error  in  the 
integrated  sectional  quantities  at  that  one  location  will  usually  destroy  the 
overall  accuracy  of  the  total  force  prediction.  Trying  to  correct  these  errors  by 
some  automatic  smoothing  process  becomes  increasingly  complex  as  the 
number  of  multiple  variables  increases  and  the  various  combinations  of 
invariant  point  motions  permutes. 

The  basic  conclusion  is  that  as  both  the  sensitivity  and  number  of  the 
base/calibration  flow  solutions  that  are  input  to  the  approximation  method 
increase,  the  accuracy  and  range  of  validity  of  the  prediction  from  the  method 
will  rapidly  reduce.  The  strategy  of  trying  to  maintain  a  constant  range  of  validity 
by  developing  increasingly  more  complex  pre-and  post-processing  procedures 
to  try  to  consistently  remove  these  inaccuracies  appears  to  be  impractical  since 
it  requires  major  verification  and  testing  and  possibly  even  redevelopment  of 
the  processing  procedures  as  both  the  number  of  design  variables  increases 
and  the  complexity  of  the  surface  flow  topology  increases.  An  alternative 
strategy  of  simply  increasing  the  numerical  accuracy  of  the  base/calibration 
solutions  by  refining  the  computational  grids  and  increasing  the  flow  field 
resolution  is  also  not  satisfactory  since  in  many  situations  it  may  not  be  feasible 
to  do  this  due  to  computer  limitations.  Furthermore,  it  has  the  negative  feature 
of  requiring  the  determination  of  higher  accuracy  flow  solutions  when  using  the 
approximation  method  than  when  not  using  it. 

in  order  to  try  to  maintain  as  wide  a  range  of  validity  of  the  approximation 
method  as  possible  as  both  the  complexity  of  the  surface  flow  topology  and 
number  of  design  variables  increases,  it  appears  that  the  most  logical  way  to 
proceed  is  to  choose  an  alternative  input  to  the  method  that  is  not  so  strongly 
sensitive  as  surface  pressures.  Flow  field  properties,  in  particular  those 
distributions  in  the  far  field,  are  a  natural  choice.  As  pointed  out  previously, 
since  many  practical  design/optimization  problems  in  both  external  and  internal 
flows  seek  improvement  of  a  figure  of  merit  represented  by  an  integrated 
quantity  such  as  drag  or  loss  coefficient,  the  information  to  determine  that 
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integrated  quantity  is  contained  as  well  in  the  far  field  flow  properties,  and  can 
be  calculated  by  employing  control  volume  integrations. 


3.5  Altemative  Formulation  of  Approximation  Method  With  Flow  Field  Rather 
Than  Surface  Properties 

We  have  posed  the  problem  of  applying  the  approximation  method  when 
employing  flow  field  information  rather  than  surface  pressure  distributions  as 
the  basic  input,  and  have  carried  out  some  preliminary  work  for  the  three- 
dimensional  problem.  At  the  full  potential  level  of  solution,  the  following  far  field 
flow  quantities  are  required:  (U,V,W)  distributions.  Since  the  flow  is  assumed 
isentropic,  pressure  and  density  distributions  can  be  determined  directly  from 
the  velocity  distributions.  At  the  Euler  and  Navier-Stokes  level  of  solution,  the 
following  far  field  quantities  are  required:  (U.V,W.RHO, ENERGY)  distributions. 
In  order  to  evaluate  the  contour  integrals  involved  in  determining  the  integrated 
forces  and  moments  acting  upon  the  configuration  inside  the  control  volume, 
the  first  step  in  the  process  is  to  determine  the  pressure  and  momentum  fluxes 
around  the  entire  outer  boundary.  For  this  purpose,  it  is  convenient  to  use  a 
contour  associated  with  the  computational  grid  at  one  or  two  mesh  points  inside 
the  outermost  mesh  boundary  location.  In  employing  the  TWING  full  potential 
code,  the  first  step  in  this  process  is  to  determine  the  physical  velocity 
components  from  the  velocity  potential.  In  the  TWING  code,  as  In  essentially  all 
current  CFD  flow  solvers,  during  the  solution  convergence  procedure  all 
computations  are  carried  out  in  a  transformed  computational  domain.  An 
illustration  of  the  overall  relationship  between  the  physical  and  computational 
domains  embodied  In  the  present  TWING  solver  is  provided  in  Figure  36. 
Physical  flow  properties  are  not  required  nor  determined  during  the 
convergence  process.  After  solution  convergence,  physical  results  are  then 
calculated  as  a  final  output  step.  Although  grid  and  potential  information  are 
available  as  output  from  the  TWING  code,  directly  differencing  the  potential  in 
physical  coordinates  will  not  retain  the  2nd  order  accuracy  maintained  in  the 
original  solution  algorithm.  In  order  to  maintain  the  2nd  order  accuracy  for  the 
physical  flow  field  distributions,  which  are  needed  for  accurately  determining 
the  total  forces,  it  is  necessary  to  difference  the  potential  in  a  fashion  consistent 
with  the  original  solution  process.  This  is  accomplished  by  employing  the  same 
grid  metrics  and  computational  molecule  surface  flux  formulations  employed  in 
the  solution  algorithm.  Although  this  information  is  embedded  in  the  flow  solver, 
extracting  the  information  directly  from  the  solver  is  not  easily  done  since 
selected  metric  arrays  are  overwritten  during  the  solution  convergence  process 
in  order  to  reduce  computer  memory  requirements.  Consequently,  we  found  it 
to  be  more  expedient  to  write  a  separate  code  to  redetermine  the  grid  metrics 
and  then  calculate  the  needed  physical  flow  field  quantities.  We  have  now 
completed  that  code.  The  computational  grid  employed  in  the  TWING  solver  is 
comprised  of  a  combination  of  two-dimensional  O-type  grids  located  in  (X.Z) 
planes  placed  at  successive  spanwise  stations  along  the  wing  from  the  wing 
root  to  beyond  the  tip,  as  illustrated  In  Figures  37  to  39  from  Ref.  10.  The  outer 
boundary  of  the  grid,  illustrated  in  perspective  in  Figure  40,  tapers  inward  at 
span  locations  along  the  wing,  and  then  becomes  cylindrical  beyond  the  wing 
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tip  out  to  the  outer  lateral  sidewall  bounda^.  Typical  results  for  physical  flow 
field  quantities  at  one  grid  point  location  inside  the  grid  outermost  boundary  are 
provided  in  Figures  41  to  43  where  distributions  of  the  velocity  components 
(U.V.W)  are  provided  at  each  of  the  30  separate  spanwise  stations.  The 
distributions  of  the  U  velocity  component  are  quite  smooth,  while  not 
surprisingly  both  the  V  an  W  components  exhibit  the  influence  of  the  wing  tip 
trailing  vortex  passing  through  the  grid  directly  downstream  of  the  wing.  The 
Influence  of  the  wing  tip  vortex  is  shown  even  more  graphically  in  Figure  44 
which  displays  density  contours  on  the  same  grid  surface. 

With  the  physical  flow  properties  in  the  region  of  the  outer  grid  boundary 
determined,  we  have  now  proceeded  to  write  the  procedure  for  evaluating  the 
pressure  and  momentum  fluxes  through  each  of  the  surfaces  defining  the 
control  volume  of  interest,  i.e.  the  quasi-cylindrical  grid  surface  enclosing  the 
wing,  the  circular  surface  boundary  at  the  symmetry  plane,  and  the 
corresponding  circular  surface  boundary  at  the  freestream  sidewall  boundary 
as  shown  in  Figure  36. 

With  the  code  complete  for  providing  second-order  accurate  physical  flow  field 
properties  from  the  TWING  solver  at  any  arbitrary  location  on  the  computational 
grid,  we  have  proceeded  to  complete  the  procedure  for  evaluating  the  pressure 
and  momentum  flux  integrals  on  a  control  surface  formed  by  the  quasi- 
cylindrical  grid  surface  enclosing  the  wing,  i.e.  a  K  s  constant  grid  shell  and  the 
the  bounding  circular  end  planes  at  the  symmetry  plane  and  the  freestream 
sidewall  boundary,  as  illustrated  in  Figure  36.  The  forces  acting  on  the  wing 
contained  within  the  control  volume  can  be  evaluated  by  determining  the 
following  two  surface  integrals: 


F  =  -J(p-Poo)dS- J| 


P(V  -  V^).(V-dS) 


(25) 


where  S  represents  the  three  surfaces  referred  to  above  and  illustrated  in 
Figure  36.  The  procedure  is  realized  numerically  by  determining  the  cell 
surface  areas  and  normal  vectors  on  the  outward-fad ng  surfaces  on  each  side 
of  the  grid  cells  contained  on  the  control  surface.  Implemetation  is  done  by 
carrying  out  the  integration  over  a  selected  K  »  Constant  radial  computational 
grid  shell  (i.e.  over  ail  I's.J's)  as  well  as  the  two  lateral  boundary  surfaces 
located  at  the  symmetry  plane  (J^i)  and  the  sidewall  boundary  (JsJMAX). 
Appropriate  control  volume  surfaces  can  range  from  the  K  =  2  grid  shell 
planned  for  use  here  down  to  the  K  »  KMAX  shell  which  conforms  to  the  wing 
and  ficticious  wing  extension  surface  (See  Fig.  36).  A  check  calculation  at  the  K 
«  KMAX  surface  should  identically  r^uce  to  the  previous  determination  of  the 
forces  by  surface  pressure  integrations  since  contributions  from  the  momentum 
flux  integral  in  Equation  (25)  should  be  zero  along  both  the  wing  and  fictidous 
wing  extension  surfaces. 
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We  have  carried  out  test  calculations  for  control  volumes  ranging  from  the  K  «  2 
to  K  =  KMAX  grid  shells.  Results  for  the  total  lift  and  drag  for  the  K  »  KMAX 
control  volume  from  our  evaluation  of  Equation  (25)  produced  identical  results 
to  those  previously  obtained  from  TWING  using  surface  pressure  integrations. 
The  momentum  flux  terms  at  the  wing  surface  were  of  the  order  (10).  and  the  net 
momentum  flux  on  the  ficticious  wing  extension  was  identically  zero.  However, 
the  results  for  the  total  forces  for  the  control  volumes  from  K  »  2  to  K  s  KMAX-1 
did  not  check  either  with  each  other  or  with  the  results  obtained  previously  from 
surface  pressure  integrations.  A  careful  check  of  the  code  for  obvious  errors  did 
not  reveal  the  source  of  the  discrepancy.  As  an  additional  check,  we  next 
examined  the  conservation  of  mass  using  the  same  control  volume  approach. 
Since  the  full  potential  equation  employed  here  is  identically  the  continuity 
equation,  conservation  of  mass  should  be  satisfied  for  any  arbitrary  control 
volume  chosen,  i.e.  from  individual  computational  grid  cells  to  the  far-field 
control  volumes  used  above  for  the  total  force  computations.  The  results  from 
the  mass  flux  integrations  over  the  same  K  s  Constant  grid  shells  plus  two 
bounding  end  planes  used  in  the  momentum  computations  produced 
erroneous  non-zero  mass  fluxes.  We  then  proceeded  to  try  to  isolate  the  source 
of  the  error.  We  decomposed  the  far-field  control  surface  into  the  sum  of 
component  spanwise  rings  formed  by  chordwise  strips  (all  I's)  between 
adjacent  spanwise  J  nodes.  We  then  examined  the  total  mass  flux  entering  and 
leaving  a  ring  control  volume  defined  by  these  individual  rings  that  were  one 
radial  K  cell  thick.  We  found  that  the  net  mass  flux  was  also  not  zero  for  these 
rings.  We  next  performed  a  check  on  the  cell  surface  area  and  cell  surface 
normal  determinations  by  imposing  a  freestream  velocity  everywhere  and  then 
computing  the  mass  flux.  A  non-zero  result  for  mass  flux  for  this  case  would 
indicate  an  error  in  either  or  both  of  the  determinations  of  the  ceil  surface  areas 
and  cell  surface  normals.  The  results  of  this  computation  was  that  zero  mass 
flux  (within  1 0  )  was  found  for:  (1 )  individual  cells,  (2)  spanwise  rings  one  K  cell 
thick,  and  (3)  the  far-field  control  volumes  used  in  the  momentum  computations. 
This  result  verified  the  computation  of  the  cell  surface  areas  and  normals.  We 
next  verified  the  accuracy  of  the  physical  (U.V.W)  velocity  distributions  being 
determined  and  employed  in  our  current  procedure.  Since  the  TWING  potential 
flow  is  isentropic,  the  velocity  distributions  can  be  used  to  determine  all  the 
subsequent  flow  properties..  We  used  our  caculated  physical  (U,V.W)  velocity 
distributions  to  determine  the  corresponding  density  distributions.  We  then 
compared  that  result  with  that  determined  by  the  TWING  solution  procedure  in 
the  computational  domain  and  which  uses  contravariant  velocity  components. 
The  density  results  from  each  of  these  separate  computations  compared 
exactly.  Finally  we  met  with  and  discussed  these  results  with  the  NASA/Ames 
Research  Center  scientists  Dr.  T.L.  Hoist  and  Mr.  S.D.  Thomas  who  are  the 
authors  of  the  TWING  code.  They  were  not  aware  of  this  discrepancy  in  the 
code,  nor  of  anyone  else  who  has  tried  to  use  far-field  distributions  from  the 
code  to  determine  total  forces.  During  development  of  the  solver,  they  had 
checked  the  predictions  of  TWING  with  a  number  of  other  3-D  solvers  and  had 
made  successful  verifications  of  pressure  distributions  and  forces  obtained  from 
the  code.  As  a  final  result  of  the  control  volume  formulation,  we  found  that  if  we 
sum  the  individual  pressure  and  momentum  flux  surface  integrations  over  all 
individual  cells  throughout  the  total  control  volume,  we  then  obtain  the  same 
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result  for  the  lift  and  drag  previously  determined  from  surface  pressure 
integrations.  This  result  implies  that  the  summing  procedure  we  are  using  is 
correct  and  that  ail  internal  fluxes  cancel  correctly. 

With  the  resources  remaining  to  the  investigation,  it  was  decided  at  this  point 
that,  rather  than  continue  to  try  to  resolve  the  discrepancies  encountered  with 
the  3-D  TWING  flow  code  in  determining  the  lift  and  drag  forces  acting  on  the 
configuration  vis-a-vis  employing  contour  integration  of  momentum  fluxes  on  an 
outer  boundary  removed  from  the  actual  geometrical  surface,  a  preliminary 
investigation  employing  flow  field  rather  than  surface  properties  with  the 
approximation  method  would  be  carried  out  at  the  axisymmetric  level.  The  flow 
solver  selected  for  use  in  this  preliminary  study  was  one  developed  by  the 
present  author  and  successfully  tested  in  an  extensive  series  of  applications 
(Refs.  30-32)  involving  highly-supercritical  nonlinear  flows  throughout  the 
transonic  regime.  The  solver  employs  a  SLOR  procedure  to  solve  the  transonic 
small-disturbance  equations  in  fully-conservative  form,  and  has  demonstrated 
its  ability  to  be  accurate  and  robust  in  determining  accurate  flow  solutions  not 
only  throughout  the  transonic  regime  but  at  subsonic  and  low  supersonic 
speeds  as  well.  The  method  has  been  configured  to  treat  free-air  as  well  as  a 
variety  of  wind-tunnel  outer  boundary  conditions  (Ref.  31).  Hence,  it  can  be 
directly  applied  to  the  present  problem  of  predicting  the  aerodynamic  forces 
acting  on  a  configuration  by  using  contour  integration  of  momentum  fluxes  over 
an  outer  control  surface.  With  this  flow  solver  we  formulated  an  initial  test  case 
for  the  approximation  method  based  on  using  flow  field  rather  than  surface 
properties  by  considering  the  axisymmetric  transonic  flow  past  a  body  of 
revolution  enclosed  by  a  solid  sidewall  boundary,  as  illustrated  in  Figure  45. 
The  control  surface  employed  for  the  contour  integration  of  the  momentum  flux 
coincided  with  the  circular  inflow  and  outflow  surfaces  at  the  upstream  and 
downstream  boundaries,  together  with  the  cylindrical  sidewall  surface.  The 
body  geometry  selected  was  chosen  to  be  a  particular  parabolic-arc  geometry 
for  which  accurate  wind  tunnel  data  are  available  from  Ref.  33.  The  particular 
parabolic-arc  body  has  a  thickness  ratio  of  1/12  together  with  a  cylindrical  sting 
attached  to  the  rear  of  the  body  at  the  85%  length  location.  For  the  calculations, 
the  sidewall  boundary  was  located  laterally  at  6  body  lengths  from  the  body 
longitudinal  axis,  and  the  upstream  and  downstream  boundary  locations  were 
set  at  3  body  lengths,  respectively,  from  the  nose  and  tail  of  the  body.  Drag 
calculations  were  then  carried  out.  both  by  integration  of  body  surface 
pressures  and  by  contour  integration  of  momentum  fluxes  over  the  control 
surface  discussed  above,  for  oncoming  Mach  numbers  ranging  from  M  =  0.90  to 
1 .20  in  steps  of  0.025.  The  agreement  between  the  two  calculations  was  very 
satisfartorily,  as  indicated  in  Figure  45.  Next,  the  approximation  method 
determination  was  reformulated  in  terms  of  the  distributions  of  the  two  velocity 
components  on  the  various  control  surfaces.  For  this  simple  problem,  this  was 
only  necessary  at  the  downstream  boundary  since  at  the  upstream  boundary, 
constant  oncoming  conditions  prevail,  while  at  the  lateral  solid  sidewall 
boundary  the  momentum  flux  is  identically  zero.  Comparative  calculations  with 
the  approximation  method  were  then  carried  out  employing  various  base  and 
calibration  solutions  to  try  to  reproduce  the  drag  curve  determined  by  the  exact 
methods.  The  particular  approximation  results  displayed  in  Figure  45  employed 
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base  and  calibration  solutions  for  M  s  {0.95.  1.025}  respectively,  but  are  typical 
of  approximation  results  obtained  for  other  choices  of  base  and  calibration 
solutions.  As  can  be  seen  from  the  figure,  the  approximation  predictions,  while 
close  to  the  exact  results,  exhibit  a  spurious  character.  We  were  unable  to 
determine  the  precise  cause  of  this  inconsistOnt  behavior,  but  suspect  that  the 
sensitivity  that  was  previously  encountered  with  the  approximation  method  in 
the  vicinity  of  shock  waves  when  using  surface  pressure  distributions  has  now 
been  transferred  to  a  corresponding  sensitivity  in  the  flow  field  velocity 
distributions  in  the  near  vicinity  of  the  body  surface.  A  corresponding  analog  for 
the  three-dimensional  problem  previously  considered  above  with  the  TWING 
code  would  be  the  sensitivity  in  the  flow  field  velocity  distributions  in  the  vicinity 
of  the  wing  trailing  vortices  that  was  noted  and  displayed  in  Figures  41  to  44. 
Successful  treatment  of  these  regions  of  high-gradient  sensitive  flow  vis-a-vis 
the  approximation  method  requires  a  reformulation  of  the  method  in  terms  of 
appropriate  invariant  points  for  whatever  flow  quantities  or  distributions  are 
selected  for  use.  This  corresponds  in  a  direct  fashion  to  the  methodology 
developed  in  this  present  study  for  the  problem  of  characterizing  various  shock 
wave  invariant  points  when  employing  surface  pressure  distributions  as  the 
basic  input  to  the  approximation  method. 
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4.  CONCLUSIONS  AND  RECOMMENDATIONS 

An  investigation  was  conducted  involving  the  preliminary  development  of  a 
three-dimensional  nonlinear  approximation  procedure  for  determining  rapid 
and  accurate  approximations  to  highly  nonlinear  flows  such  as  typically  occur  in 
turbomachinery  applications.  The  overall  objective  was  to  develop  and 
demonstrate  the  feasibility  of  such  approximation  methods  to  successfully  work 
in  a  general  3-D  turbomachinery  design  or  parametric  optimization  environment 
to  substantially  reduce  the  overall  computational  requirements  necessary  for 
optimum  design.  The  approximation  procedures  employ  unit  perturbations, 
determined  from  two  or  more  nonlinear  "base”  solutions  which  differ  from  one 
another  by  a  nominal  change  in  some  geometry  or  flow  parameter,  to  predict  a 
family  of  related  nonlinear  solutions.  The  solutions  can  be  either  continuous  or 
highly  discontinuous.  It  is  conceived  that  these  methods  would  be  coupled  with 
high  run-time  three-dimensional  computational  flow  solvers,  and  would  be  used 
in  conjunction  with  these  solvers  in  applications  where  large  numbers  of  related 
solutions  are  needed.  The  computational  time  saving  would  be  accomplished 
by  employing  these  rapid  approximation  methods  to  decrease  the  actual 
number  of  expensive  3-D  flow  solutions  needed  in  any  optimization  or  design 
study  to  a  minimum. 

The  work  undertaken  here  relates  to  the  initial  development  and  extension  of 
these  methods  and  concepts  to  problems  characteristic  of  three-dimensional 
turbomachinery  optimization  design.  The  specific  objectives  of  this  study  were: 
theoretical  development  of  the  three-dimensional  approximation  procedure  in  a 
form  suitable  for  predicting  surface  properties  on  three-dimensional 
turbomachinery  blades  for  highly  sensitive  supercritical  transonic  flows 
involving  multiple  shock  waves,  combination  of  the  approximation  method  with 
a  nonlinear  three-dimensional  transonic  flow  solver,  coupling  of  the  combined 
approximation  method  and  nonlinear  flow  solver  with  an  optimization  design 
procedure,  and  finally  testing  of  the  complete  approximation  method/3-D  flow 
solver/optimization  code  on  problems  relevant  to  the  three-dimensional 
turbomachinery  design  optimization  environment. 

Application  of  the  approximation  method  to  three-dimensional  transonic  flows 
with  multiple  shock  waves  required  the  development  of  a  multiple  invariant 
point  characterization  procedure  for  identifying  and  tracking  through  parameter 
solution  space  the  characteristic  topology  features  associated  with  shock  waves 
appearing  in  the  surface  pressure  distributions  used  as  input  to  the 
approximation  method.  It  was  found  from  numerical  experimentation  that  a  two- 
point  characterization  of  shock  topology  in  surface  pressure  distributions  was 
adequate  for  the  transonic  flows  of  interest  in  this  study.  The  characterization 
procedure  was  extended  to  include  additional  invariant  points  required  for 
topology  tracking  of  other  high-gradient  features  characteristic  of  the  general 
class  of  three-dimensional  transonic  surface  pressure  distributions  being 
considered.  These  include  stagnation  point  and  minimum  suction  pressure  and 
associated  recovery  points  near  the  nose  on  both  upper  and  lower  surfaces. 
Finally,  a  global  invariant  point  preprocessing  procedure  was  developed  to 
ensure  a  continuous  variation  in  the  spanwise  coordinate  direction  of  all  the 
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selected  invariant  points.  The  approximation  method  was  then  tested  on  a 
number  of  three-dimensional  example  problems  involving  strongly  supercritical 
transonic  flows.  The  three-dimensional  transonic  full  potential  flow  solver 
TWING  was  employed  to  provide  the  flow  solutions.  Comparative  results  of  the 
approximation  method  predictions  with  the  exact  nonlinear  solutions  indicated 
that  the  approximation  method  is  able  to  provide  good  nonlinear  predictions 
over  relatively  large  parametric  ranges  for  three-dimensional  problems 
involving  one  to  three  simultaneous  parameter  variations. 

Existing  aerodynamic  design  optimization  methods  and  applications  thereof 
were  reviewed  and,  as  a  result,  the  quasi-Newton  optimization  driver  QNMDIF 
was  selected  as  superior  for  optimization  studies  of  the  general  class  of  interest 
to  this  study.  That  optimization  method  was  then  coupled  with  the  three- 
dimensional  approximation  method  and  the  three-dimensional  TWING  flow  field 
solver.  Using  past  studies  with  the  QNMDIF  optimization  driver  as  a  guide,  a  3- 
D  design  optimization  problem  was  then  defined  for  testbed  case  studies 
involving  the  combined  code.  The  testbed  optimization  problem  selected 
employed  as  design  variables  the  vertical  displacement  of  a  number  of  pre¬ 
assigned  points  located  in  a  certain  region  on  the  geometric  surface  of  a  3-D 
isolated  blade  which  was  to  be  modified,  and  incorporated  a  spline-fitting 
procedure  to  redefine  the  surface  geometry  as  the  pre-assigned  points  displace 
vertically  from  their  baseline  locations  in  response  to  the  optimization  process. 
The  drag/lift  ratio  was  selected  for  the  objective  function.  The  testbed  problem 
was  first  solved  using  only  full  nonlinear  3-D  TWING  flow  solutions  in  order  to 
obtain  the  exact  optimum  design  point  toward  which  the  approximation  method 
could  be  tested.  Application  of  the  approximation  method  was  then  made  to  the 
testbed  problem  to  examine  the  capability  of  the  method  to  produce  accurate 
results  in  a  typical  design  environment,  and  also  to  evalutate  its  potential  for 
computational  savings.  Sensitivity  studies  were  performed  to  examine  the 
accura^  dependence  of  the  approximation  method  on  the  choice  of  the  initial 
calibration  solution  matrix. 

Comparisons  of  the  three-dimensional  approximation  method,  configured  with  a 
multiple  invariant  point  characterization  for  surface  shock  waves,  stagnation 
point,  and  upper  and  lower  surface  suction  pressure  and  recovery  points,  with 
the  corresponding  exact  nonlinear  solutions  for  surface  pressure  distributions 
Indicated  good  accuracy  and  range  of  validity  of  the  approximation  method  for 
those  situations  where  the  basic  flow  topology  did  not  significantly  change  over 
the  parameter  range  studied.  For  the  sensitive  transonic  flows  considered  in 
the  testbed  optimization  case  study,  however,  it  was  found  necessary  to 
develop  an  additional  postprocessing  procedure  to  relocate  all  the  invariant 
points  in  order  to  provide  a  globally  continuous  variation  in  both  the  spanwise 
invariant  point  locations  as  well  as  their  associated  surface  pressure  values. 
With  this  enhancement  of  the  approximation  method,  the  method  was  applied  to 
predict  both  surface  pressure  distributions  and  spanwise  sectional  lift  and  drag 
distributions  within  the  9  design  variable  solution  space  encompassed  by  the 
testbed  problem.  The  comparative  results  for  the  surface  pressure  distributions 
provided  reasonable  overall  agreement  with  the  exact  results,  but  exhibited 
discrepancies  related  to  the  rapid  topology  changes  occurring  in  the  spanwise 
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direction  in  these  flows.  However,  the  spanwise  sectional  lift  and  drag 
distribution  comparisons  revealed  a  major  shortcoming  in  the  approximation 
method.  That  is,  while  the  sophisticated  invariant  point  relocation  procedures 
developed  here  insure  continuous  variations  of  invariant  point  locations  and 
surface  pressure  values  at  those  locations,  this  by  itself  does  not  insure  a 
continuous  variation  in  the  sectional  integrated  quantities  obtained  from  those 
pressure  distributions.  The  reason  for  this  is  again  due  to  the  tendency  for 
extremely  rapid  local  topology  changes  to  occur  in  the  basic  surface  pressure 
distribution  for  these  strongly  nonlinear  flows.  These  rapid  topology  changes  in 
surface  pressure  resulted  in  a  small  displacement  error  in  the  spanwise 
invariant  point  location  distribution  causing  the  strength  and  location  of  the 
shock  at  one  or  more  spanwise  locations  to  be  sufficiently  in  error  that  a 
significant  error  in  the  spanwise  variation  of  the  integrated  sectional  quantities 
resulted.  To  investigate  this  problem  further,  a  detailed  evaluation  of  the 
performance  of  the  approximation  method  was  carried  out  involving  problems 
having  one  to  nine  parametric  variables.  The  following  four  major  steps  in  the 
approximation  procedure  were  carefully  reviewed:  (1)  the  preprocessing 
procedure  for  identifying  and  relocating  the  invariant  point  locations  in  the  base 
and  calibration  solutions,  (2)  the  procedure  for  determining  the  approximation- 
predicted  solution  based  upon  the  preprocessed  base/calibration  solutions,  (3) 
the  post-processing  procedure  for  relocating  the  predicted  invariant  point 
locations  in  the  spanwise  direction,  (4)  the  integration  of  the  post-processed 
approximation  predicted  solutions  to  determine  integrated  spanwise  sectional 
quantities.  As  a  result  of  the  review,  the  multiple  sources  causing  the  observed 
limitation  of  the  approximation  method  in  three-dimensional  studies  were 
identified.  The  primary  factors  involve  a  combination  of  inherent  numerical 
inaccuracies  in  the  three-dimensional  base/calibration  solutions  near  the 
invariant  point  locations  associated  with  the  foot  of  the  shock,  together  with  a 
geometrical  compounding  of  these  errors  when  the  approximation  predicted 
invariant  point  locations  are  dependent  upon  the  simultaneous  variation  of  a 
large  number  of  multiple  parameters.  The  fundamental  action  of  the 
approximation  method  is  to  move  the  invariant  points  linearly  in  multiply- 
strained  coordinate  space  with  respect  to  each  perturbation  parameter.  The 
results  predicted  from  such  a  formulation  have  been  found  to  be  extremely 
accurate  in  two-dimensional  situations  where  the  shock  structure  is  both 
relatively  simple  and  well  resolved  by  the  numerical  solution  procedure  so  that 
the  associated  invariant  points  are  also  well  defined.  For  three-dimensional 
situations,  however,  the  shock  structure  on  or  near  the  surface  of  a  3-D 
configuration  is  not  so  clearly  defined  by  the  flow  solvers.  Consequently,  the 
definition  of  the  invariant  point  locations  in  the  base/calibration  solutions  is 
correspondingly  lacking.  Attempting  to  correct  these  locations  by  employing 
global  smoothing  criteria  does  improve  the  situation  at  most  locations  but 
inevitably  increases  the  error  at  other  locations.  The  basic  conclusion  is  that  as 
both  the  sensitivity  and  number  of  the  base/calibration  flow  solutions  that  are 
input  to  the  approximation  method  increase,  the  accuracy  and  range  of  validity 
of  the  prediction  from  the  method  will  rapidly  reduce.  The  strategy  of  trying  to 
maintain  a  constant  range  of  validity  by  developing  increasingly  more  complex 
pre-and  post-processing  procedures  to  try  to  consistently  remove  these  inherent 
inaccuracies  appears  to  bo  impractical  since  it  requires  major  verification  and 
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testing  and  possibly  even  redevelopment  of  the  processing  procedures  as  both 
the  number  of  design  variables  increases  and  the  complexity  of  the  surface  flow 
topology  increases. 

It  appears  that  the  most  logical  way  to  proceed  is  to  employ  an  alternative  input 
to  the  approximation  method  that  is  not  so  strongly  sensitive  and  subject  to 
rapid  toplogy  changes  as  surface  pressures.  In  order  to  test  this  hypothesis,  the 
approximation  method  was  reformulated  in  terms  of  using  flow  field  velocity 
distributions  at  or  near  the  outer  grid  boundary,  and  an  initial  preliminary 
investigation  was  carried  out  employing  the  TWING  flow  field  solver. 
Discrepancies  encountered  in  the  mass  and  momentum  flux  determinations 
over  a  control  surface  located  near  the  outer  boundary  of  the  TWING 
computational  mesh  prevented  a  definitive  conclusion  regarding  the  accuracy 
of  the  three-dimensional  approximation  method  based  on  input  flow  field 
properties.  Initial  results  from  a  corresponding  axisymmetric  flow  case  study 
provided  promising  results  for  wave  drag  prediction,  and  pointed  to  the  need  for 
the  appropriate  characteristic  invariant  point  development  associated  with  the 
particular  flow  field  properties  employed. 

Finally,  the  results  of  the  present  investigation  provide  an  important  guideline 
for  future  development.  We  believe  that  the  demonstrated  potential  of  the 
approximation  method  from  both  the  present  investigation  and  past  studies,  in 
particular  for  design  optimization  studies  involving  highly  nonlinear  two- 
dimensional  transonic  internal  and  external  flows  where  Its  ability  to  reduce  the 
computational  work  by  an  order  of  magnitude  with  no  degradation  in  accuracy 
was  clearly  shown,  warrants  further  development  of  the  method  for  the  three- 
dimensional  turbomachinery  design  problem.  Furthermore,  since  the 
approximation  method  does  not  depend  upon  any  particular  flow  analysis  code, 
obsolescence  of  the  methodology  developed  due  to  future  analysis  code 
improvement  will  not  occur.  Flow  field  property  formulations  of  the 
approximation  method  should  be  pursued  as  being  most  suitable  for  complex 
multiple  parameter  design  problems. 
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Perturbation  for  Perturbation  for 

calibration  solution  calibration  solution 

in  physical  coordinates  in  strained  coordinates 


(a)  Single  shock 


(b)  Multiple  shock  and  high-gradient  regions 


Figure  1.  Illustration  of  perturbation  solution  for  calibration  solution  in 
physical  and  strained  coordinates 
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Figure  2.  Comparison  of  perturbation  and  exact  nonlinear  surface  pressure 
distributions  for  a  thickness-ratio  perturbation  of  an  isolat^  NACA 

OOXX  blade  profile  at  M«,  =  0.820  and  a  =  0  ®  for  solution 
interpolation 
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Comparison  of  approximation  predicted  and  exact  nonlinear 
surface  pressure  distributions  for  a  thickness-ratio  perturbation  of  a 
NACA  OOXX  blade  profile  at  M«  «  0.820  and  a  «  0  °  for  illustrating 
both  solution  interpolation  and  extrapolation 


at, -1.00 

OCg  =  2.00® 


a  »  2.50 


a  =  1.50 


■  —  Evaei 

•  •••  f^apraK.  na4hod 

Saaa 

_ ...  Callbra4laA 


a.  0.50 


Caaei 

Aa^raa.  Mathatf 
Baaa 

Cal Ibraf ian 


Comparison  of  approximation  predicted  and  exact  nonlinear  upper 
and  lower  surface  pressure  distributions  for  an  angle-of-attack 
perturbation  of  a  NACA  0012  blade  profile  at  M«  «  0.750  illustrating 
both  solution  interpolation  and  extrapolation 
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Figure  6.  Comparison  of  approximation  predicted  and  exact  nonlinear 
surface  pressure  distributions  for  a  simultaneous  perturbation  of  a 
oncoming  Mach  number  and  thickness  ratio  of  a  symmetric  NACA 
OOXX  blade  profile  illustrating  both  solution  interpolation  and 
extrapolation 


Illustration  of  the  shock  wave  sharpening  achieved  by  the 
preprocessing  procedure  involving  the  invariant  points  associated 
with  the  surface  shock  waves  on  the  3-D  supercritical  surface 
pressure  distribution  shown  in  Figure  7 


Illustration  of  the  smooth  global  spanwise  correction  achieved  by 
the  preprocessing  procedure  acting  on  all  of  the  invariant  points 
associated  with  the  3-D  supercritical  surface  pressure  distribution 
shown  in  Figure  7 


Figure  1 1 .  Comparison  of  3-D  approximation  predicted  and  exact  nonlinear 
surface  pressure  distributions  at  the  root  chord  station  (y/s  =  0.0)  for 
a  thickness  ratio  perturbation  of  the  ONERA  M6  blade  profile  shown 
in  Figure  10  illustrating  both  solution  interpolation  and 
extrapolation 
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Figure  12.  Comparison  of  3-D  approximation  predicted  and  exact  nonlinear 
surface  pressure  distributions  at  the  outermost  spanwise  station 
(y/s  s  0.969)  for  a  thickness  ratio  perturbation  of  the  ONERA  M6 
blade  profile  shown  in  Figure  10  illustrating  both  solution 
interpolation  and  extrapolation 
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Figure  13.  Illustration  of  the  extreme  surface  pressure  topology  change 
involved  in  the  oncoming  Mach  number  perturbation  test  problem; 
surface  pressure  distributions  at  the  beginning,  middle,  and  end  of 
the  parameter  range  M  =  {0.880.  0.800,  0.720} 
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Figure  14.  Comparison  of  3-D  approximation  predicted  and  exact  nonlinear 
surface  pressure  distributions  at  the  outermost  spanwise  station 
(y/s  -  0.969)  of  the  ONERA  M6  blade  profile  shown  in  Rgure  10  for 
an  oncoming  Mach  number  perturbation  illustrating  both  solution 
interpolation  and  extrapolation 


0.0  0.1  0.2  0.3  0.<  0.5  0.6 

CHORD  LENGTH 


0.7 


0.8  0.9  1.0 


Location  ot  the  fixed  (+)  and  Movable  (*)  apllna-aupport  polnta  on  the  HACA  6AjA212  airfoil 


Figure  1 5.  Outline  of  overall  optimization  design  test  problem  to  be  ernployed 
as  the  benchmark  case  study  for  testing  the  3-D  approximation 
method  In  a  design  environment;  illustration  of  the  isolated  blade 
planform  geometry  and  location  of  design  variable  points 
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Figure  16.  Illustration  of  the  3-D  surface  pressure  distribution  for  the  original 
(nonoptimized)  configuration  for  the  benchmark  optimization  case 
study 
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Figure  17.  Illustration  of  the  3-D  surface  pressure  distribution  for  the 
optimization-modified  configuration  after  the  7th  optimization  cycle 
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Figure  18.  Illustration  of  the  3-D  surface  pressure  distribution  for  the 
optimization-modified  configuration  after  the  8th  optimization  cycle 
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25  Span  S+a+ion  from  ROOT  +o  TIP 


Figure  19.  Illustration  of  the  3-D  surface  pressure  distribution  for  the 
optimization-modified  configuration  after  the  9th  optimization  cycle 
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Figure  20.  Illustration  of  the  3-D  surface  pressure  distribution  for  the 
optimization-modified  configuration  after  the  10th  optimization  cycle 
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Figure  21 .  Convergence  history  of  the  optimization  problem  described  in 
Figures  15  to  20;  plot  of  objective  function  at  each  evaluation  by  3- 
D  flow  solver  as  called  by  optimization  driver  from  problem  start 
through  the  1 0th  optimization  cycle 
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CfiSE  0-3:  3D  ONM/TWING  PERTURBRTION 
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Figure  22.  Comparison  of  3-0  approximation  predicted  and  exact  nonlinear 
surface  pressure  distributions  at  the  root  chord  station  (y/s  =  0.0)  for 
the  particular  values  of  the  9  design  variables  at  the  end  of  the  7th 
optimization  cycle;  complete  exact  nonlinear  3-D  surface  pressure 
distribution  shown  in  Figure  17 
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Figure  23.  Comparison  of  3-D  approximation  predicted  and  exact  nonlinear 
surface  pressure  distributions  at  the  root  chord  station  (y/s  =  0.0)  for 
the  particular  values  of  the  9  design  variables  at  the  end  of  the  8th 
optimization  cycle;  complete  exact  nonlinear  3-D  surface  pressure 
distribution  shown  In  Figure  18 
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Figure  24.  Comparison  of  3-D  approximation  predicted  and  exact  nonlinear 
surface  pressure  distributions  at  the  root  chord  station  (y/s  =  0.0)  for 
the  particular  values  of  the  9  design  variables  at  the  end  of  the  9th 
optimization  cycle;  complete  exact  nonlinear  3-D  surface  pressure 
distribution  shown  in  Figure  19 
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Figure  25.  Comparison  of  3-D  approximation  predicted  and  exact  nonlinear 
surface  pressure  distributions  at  the  root  chord  station  (y/s  =  0.0)  for 
the  particular  values  of  the  9  design  variables  at  the  end  of  the  10th 
optimization  cycle;  complete  exact  nonlinear  3-D  surface  pressure 
distribution  shown  in  Figure  20 
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Figure  26.  Illustration  of  the  spanwise  sectional  lift,  drag,  and  drag/lift  ratio  at 
each  of  the  25  spanwise  locations  along  the  3-D  profile  for  the 
surface  pressure  distribution  shown  as  obtained  at  the  end  of  the 
7th  optimization  cycle  for  the  benchmark  case  study 


Figure  27.  Illustration  of  the  spanwise  sectional  lift,  drag,  and  drag/lift  ratio  at 
each  of  the  25  spanwise  locations  along  the  3-D  profile  for  the 
surface  pressure  distribution  shown  as  obtained  at  the  end  of  the 
8th  optimization  cycle  for  the  benchmark  case  study 
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Figure  29.  Illustration  of  the  spanwise  sectional  lift,  drag,  and  drag/lift  ratio  at 
each  of  the  25  spanwise  locations  along  the  3'D  profile  for  the 
surface  pressure  distribution  shown  as  obtained  at  the  end  of  the 
1 0th  optimization  cycle  for  the  benchmark  case  study 


Figure  30.  Illustration  of  improvement  in  spanwise  invariant  point  location  and 
invariant  point  Cp  values  provided  by  invariant  point  relocation 

procedure  incorporated  in  3-0  approximation  method  for  baseline 
3-D  surface  pressure  distribution  for  benchmark  case  study 


Figure  31 .  Illustration  of  improvement  in  spanwise  invariant  point  location  and 
invariant  point  Cp  values  provided  by  invariant  point  relocation 

procedure  incorporated  in  3-0  approximation  method  for  3*D 
surface  pressure  distribution  calibration  associated  with  calibration 
solution  for  2nd  design  variable  for  benchmark  case  study 


Figure  32.  (lustration  of  improvement  in  spanwise  invariant  point  location  and 
invariant  point  Cp  values  provided  by  invariant  point  relocation 
procedure  incorporated  in  3-0  approximation  method  for  3-D 
surface  pressure  distribution  calibration  associated  with  calibration 
solution  for  5th  design  variable  for  benchmark  case  study 
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Figure  33.  Comparison  of  3-D  approximation  predicted  and  exact  nonlinear 
chordwise  surface  pressure  distributions  at  ail  25  spanwise  stations 
for  the  pressure  distribution  associated  with  the  7th  optimization 
cycle  design  result  for  the  benchmark  case  study 
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Figure  33.  -  Continued 
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Figure  33.  -  Continued 
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Figure  34.  Comparison  of  3-0  approximation  predicted  and  exact  nonlinear 
spanwise  distributions  of  the  sectional  lift,  drag,  and  drag/lift  ratio 
associated  with  the  7th  optimization  cycle  design  result  for  the 
benchmark  case  study 


Figure  35.  (lustration  of  the  approximation  predicted  spanwise  distributions  of 
the  sectional  lift,  drag,  and  drag/lift  ratio  based  on  post'processed 
invariant  p}oint  surface  pressure  distributions 


Figure  36. 
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llustration  of  the  physical  and  computational  computational 
domains  embodied  in  the  3-D  TWING  flow  field  solver 
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Figure  37.  Ilustration  of  the  two-dimensional  (X.Z)  0-type  grid  at  the  first 
spanwise  station  (symmetry  plane)  along  wing 
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Figure  38.  Hustration  of  the  two-dimensional  (X,Z)  0-type  grid  at  the  first 
spanwise  station  beyond  the  wing  tip 
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Figure  39.  Planform  view  of  TWING  grid  in  Z  =  0  plane  illustrating  grid  tapering 
and  clustering  near  leading  and  trailing  edges  and  ficticious  wing 
extension 
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Figure  40.  Perspective  view  of  outer  boundary  of  3-D  TWING  grid;  K  «  1 
computational  grid  shell  in  physical  coordinates 
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41 .  Distribution  of  U  velocity  component  at  the  30  spanwise  locations 
from  wing  root  to  sidewall  boundary  on  the  K  »  2  grid  shell 
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Figure  43.  Distribution  of  W  velocity  component  at  the  30  spanwise  locations 
from  wing  root  to  sidewall  boundary  on  the  K  »  2  grid  shell 
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Figure  44,  Density  contours  on  the  K  =  2  computational  grid  shell  illustrating 
the  effect  of  the  wing  tip  trailing  vortex 
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Figure  45.  illustration  of  test  case  formulation  for  approximation  method  based 
on  using  flow  field  rather  than  surface  properties;  axisymmetric 
transonic  flow  past  a  body  of  revolution  with  a  solid  wall  outer 
boundary ;  comparison  of  wave  drag  predictions 
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